Chromatin Interaction Responds to Breast Muscle Development and Intramuscular Fat Deposition Between Chinese Indigenous Chicken and Fast-Growing Broiler

Skeletal muscle development and intramuscular fat (IMF) content, which positively contribute to meat production and quality, are regulated by precisely orchestrated processes. However, changes in three-dimensional chromatin structure and interaction, a newly emerged mediator of gene expression, during the skeletal muscle development and IMF deposition have remained unclear. In the present study, we analyzed the differences in muscle development and IMF content between one-day-old commercial Arbor Acres broiler (AA) and Chinese indigenous Lushi blue-shelled-egg chicken (LS) and performed Hi-C analysis on their breast muscles. Our results indicated that significantly higher IMF content, however remarkably lower muscle fiber diameter was detected in breast muscle of LS chicken compared to that of AA broiler. The chromatin intra-interaction was prior to inter-interaction in both AA and LS chicken, and chromatin inter-interaction was heavily focused on the small and gene-rich chromosomes. For genomic compartmentalization, no significant difference in the number of B type compartments was found, but AA had more A type compartments versus LS. The A/B compartment switching of AA versus LS showed more A to B switching than B to A switching. There were no significant differences in the average sizes and distributions of topologically associating domains (TAD). Additionally, approximately 50% of TAD boundaries were overlapping. The reforming and disappearing events of TAD boundaries were identified between AA and LS chicken breast muscles. Among these, the HMGCR gene was located in the TAD-boundary regions in AA broilers, but in TAD-interior regions in LS chickens, and the IGF2BP3 gene was located in the AA-unique TAD boundaries. Both HMGCR and IGF2BP3 genes exhibited increased mRNA expression in one-day-old AA broiler breast muscles. It was demonstrated that the IGF2BP3 and HMGCR genes regulated by TAD boundary sliding were potential biomarkers for chicken breast muscle development and IMF deposition. Our data not only provide a valuable understanding of higher-order chromatin dynamics during muscle development and lipid accumulation but also reveal new insights into the regulatory mechanisms of muscle development and IMF deposition in chicken.

Skeletal muscle development and intramuscular fat (IMF) content, which positively contribute to meat production and quality, are regulated by precisely orchestrated processes. However, changes in three-dimensional chromatin structure and interaction, a newly emerged mediator of gene expression, during the skeletal muscle development and IMF deposition have remained unclear. In the present study, we analyzed the differences in muscle development and IMF content between one-dayold commercial Arbor Acres broiler (AA) and Chinese indigenous Lushi blue-shelled-egg chicken (LS) and performed Hi-C analysis on their breast muscles. Our results indicated that significantly higher IMF content, however remarkably lower muscle fiber diameter was detected in breast muscle of LS chicken compared to that of AA broiler. The chromatin intra-interaction was prior to inter-interaction in both AA and LS chicken, and chromatin inter-interaction was heavily focused on the small and gene-rich chromosomes. For genomic compartmentalization, no significant difference in the number of B type compartments was found, but AA had more A type compartments versus LS. The A/B compartment switching of AA versus LS showed more A to B switching than B to A switching. There were no significant differences in the average sizes and distributions of topologically associating domains (TAD). Additionally, approximately 50% of TAD boundaries were overlapping. The reforming and disappearing events of TAD boundaries were identified between AA and LS chicken breast muscles. Among these, the HMGCR gene was located in the TAD-boundary regions in AA broilers, but in TADinterior regions in LS chickens, and the IGF2BP3 gene was located in the AA-unique TAD boundaries. Both HMGCR and IGF2BP3 genes exhibited increased mRNA expression in one-day-old AA broiler breast muscles. It was demonstrated that the IGF2BP3 and HMGCR genes regulated by TAD boundary sliding were potential biomarkers for chicken breast muscle development and IMF deposition. Our data not only provide a valuable understanding of higher-order chromatin dynamics during muscle development

INTRODUCTION
Chicken is often regarded as one of the most desirable meats worldwide due to its proteins, polyunsaturated fatty acids (PUFA), calcium, phosphorus, and iron content, as well as its potential medicinal function and low cost (Soriano-Santos, 2010). In recent decades, the genetic selection for body weight, growth rate, and feed conversion rate of broiler chickens has contributed to a prominent increase in production efficiency, such as meat yield; however, there has also been a parallel decrease in meat quality. The synchronized improvement of both meat yield and quality has become a hot topic in molecular breeding in chicken. Skeletal muscle constitutes the largest proportion and most valuable component of meat yield, and its development and composition are of extreme importance in meat production efficiency and quality, particularly breast muscle in poultry.
Skeletal muscle is mainly composed of muscle fibers derived from the first formation of primary muscle fibers at 3.5 embryos old and the subsequent formation of secondary muscle fibers at 15 embryos old until hatching during the embryonic period in chicken (Bandman and Rosser, 2000;Picard et al., 2002). Myoblasts undergo first proliferation, followed by the fusion to form myotubes, then final differentiation into muscle fibers. Muscle growth and development go through two stages, an increase in the number of proliferative myofibers (termed hyperplasia) from the embryonic stage to early growth stage, and an increase in the size of myofibers (termed hypertrophy) during late growth stages. It has been demonstrated that the total number and morphological structure of myofibers are fixed before hatching or birth in poultry (Smith, 1963). Muscle development after birth mainly profits from growing muscle fiber diameter and length, which involves complex regulatory networks and signaling pathways. Compared to slow-growing native chicken, fast-growing commercial broilers showed a significant increase in diameter, cross-sectional area, and density of breast muscle fibers at hatching. It has been proposed that heavier chicks could have enhanced satellite cell proliferation and earlier differentiation. The increased activity of muscle satellite cells in the late stage of the embryo could make breast muscle development faster (Velleman et al., 2000;Sklan et al., 2003).
It was reported that myocytes and adipocytes could affect each other during development, and meat-producing animals that have increasing muscularity displayed a reduced intramuscular fat (IMF) development (Hocquette et al., 2010). IMF refers to the lipid between muscle fibers, mainly distributed in epimysium, perimysium, and endomysium. IMF content serves as a deeply predominant factor positively underlying meat quality, including flavor, tenderness, and juiciness, since IMF is particularly rich in phospholipids, an important precursor substance carrying unsaturated fatty acid. IMF content depends upon the hyperplasia (number) and hypertrophy (size) of adipocytes in muscle. Generally, the number of adipocytes, which are mainly determined by genetic factors, is stationary during embryonic and early development stages. The size of adipocyte refers to the ability of lipid droplet accumulation, which is determined by lipid hydrolysis and traffic of both endogenous exogenous triglyceride in the blood, uptake across cytomembrane, and metabolism including fatty acids esterification, oxidation, and lipolysis in the cytoplasm (Qiu et al., 2017). IMF deposition is mainly regulated by genetic and nutritional factors with medium to high heritability. Increasing research on the difference in meat quality among chicken breeds has emphasized the comparison of fastgrowing commercial broilers and slow-growing native chicken breeds, indicating that the latter has superiority and priority in meat quality and palatability (Leenstra, 1986;Chen et al., 2005;Li et al., 2006;Hocquette et al., 2010;Chen et al., 2016). It has been established that the IMF level at the beginning of the growth period was a likely key determinant of the final IMF level in cattle, indicating such significance of IMF accumulation in the early period. The chicken IMF content showed a significant increase from day 56 (fast-growth age), day 98 (marketing age), and day 140 (first egg age), but reached a peak at hatching in which the protein expression profile of breast muscles was especially varied from the other ages (Liu et al., 2016). It was similar to duck, which manifests a maximum lipid level and the relative area occupied by adipocytes on cross-sections of breast muscle (Chartrin et al., 2007). Another case in point also showed that the major period of IMF deposition in chicken was from embryo day 17 to day 1 (Liu R. et al., 2017). This evidence supports the idea that more attention should be paid to IMF deposition at the early period, especially the hatching.
It has been confirmed that no genetic conflict occurs between meat quality and meat production, implying that selection is valuable to synchronously improve them in broilers (Le Bihan-Duval et al., 2008). Both skeletal muscle development and IMF deposition are sophisticated and precisely orchestrated processes mediated by a network of regulatory factors. Although several studies have highlighted the roles of coding and noncoding genes in skeletal muscle development and IMF deposition, their precise molecular mechanism has largely remained elusive, especially three-dimensional chromatin structure and interaction, which have been gradually characterized as critical new regulators and have been recognized as an emerging mediator for triggering gene expression (de Laat and Duboule, 2013;Dekker and Mirny, 2016).
High-through chromosome conformation capture (Hi-C) technique has recently been developed to study the genomewide spatial chromatin organization and obtain genome-wide and high-resolution chromatin intra-(cis) and inter-(trans) interaction to elucidate gene regulation mechanisms, consequently phenotype formation (Belton et al., 2012;Stadhouders et al., 2012;Liu et al., 2019). Chromatin is characterized with high conformation by hierarchical structural units organized into following chromosome territory (CT), compartments including transcriptionally active gene-rich compartment A (open euchromatin) and transcriptionally inactive gene-desert compartment B (closed heterochromatin), regions of more localized interaction known as topologically associating domain (TAD) as well as chromatin loop in the nucleus (Gibcus and Dekker, 2013). It has been largely acknowledged that compartment B has a higher interaction intensity than compartment A, and A/B compartment switch (A to B and B to A) directly aids in the activation or inhibition of gene expression. TADs represent local contact-enriched areas within the compartment, which correspond to relatively isolated packing units of gene expression regulation with the spatial proximity of genes and their regulatory elements located in a long linear region (Dixon et al., 2012). Between TADs, there are boundaries, called TAD boundaries. TAD boundaries maintain TAD independence to characterize by self-association and insulation, which respectively promotes the gene promoter interaction with regulatory elements within TAD or avoids that across TAD borders, and play a crucial role in stabilizing genome structure and regulating gene expression (de Laat and Duboule, 2013;Ibrahim and Mundlos, 2020). However, the disappearing and remodeling of TAD boundaries may lead to abnormal intergenic interactions and altered gene expression levels (Dixon et al., 2015). Chromatin loops facilitate long-range interactions of the promoter with enhancer and silencer. They can affect gene regulation to a certain extent while dynamic formation and the disappearance of the loop structure occur (Stadhouders et al., 2012).
Lushi blue-shelled-egg chicken (LS) is a Chinese indigenous breed with the strengths of crude feed tolerance, strong adaptability, high disease resistance, as well as favorable meat quality. However, this breed has a low growth rate, feed conversion ratio, and meat production. Arbor Acres broiler (AA) is a specialized fast-growing commercial broiler, characterized as having fast weight gain, a high feed conversion ratio, and massive meat production, especially breast muscle.
To determine the regulation mechanisms underlying genome organization and genome-wide chromatin interactions during skeletal muscle development and IMF deposition, we comprehensively determined the chromatin interaction, A/B compartment switch, TAD distribution, and potential genes and signaling pathways, and unveiled gene regulation mediated by chromatin interaction at the chromosome scale using Hi-C in the breast muscles from commercial broiler AA and Chinese native Lushi chicken. Our data not only reveal a connection between three-dimensional chromatin interaction with muscle development and IMF deposition but also provide a valuable genomic resource for new insights into understanding the molecular mechanisms of skeletal muscle development and IMF deposition in chicken.

Ethics Statement
All birds used in this study were obtained from the Animal Center of Henan Agricultural University. The procedures of all animal experiments were approved by the Institutional Animal Care and Use Committee of Henan Agricultural University (Permit Number: 11-0085). All efforts were made to humanely slaughter the birds and consider animal welfare.

Experimental Animals and Sample Preparation
All the experimental chickens were raised in the same feeding environment with free access to water and feed. Breast muscles tissues were collected from 10-embryonic-old (E10), 14embryonic-old (E14), 18-embryonic-old (E18), 1-day-old (1D), 3-week-old (3W), and 5-week-old (5W) AA and LS chicken (n 6), respectively. Subsequently, they were immediately frozen in liquid nitrogen and stored at −80°C. A total of 30 5-week-old male Lushi and AA chickens (15 individuals for each breed) were selected to measure body weight, shank length, shank circumference, sternum length, body oblique length, as well as depth, width, and angle of breast muscle. Then, blood was collected from the wing vein into a 5 ml coagulation tube, underwent an incubation for 2 h, and centrifuged at 3,000 r/ min for 10 min at room temperature for serum collection. According to the manufacturer's specifications (Huili, Changchun, China), the serum biochemical indexes including alanine transaminase (ALT), aspartate transaminase (AST), glucose (GLU), lactate dehydrogenase (LDH), triglyceride (TG), total cholesterol (T-CHO), high density lipoprotein (HDL-c), low density lipoprotein (LDL-c) and very lowdensity lipoprotein (VLDL-c) were tested. After being humanely slaughtered, the breast muscle was harvested and infiltrated in 4% paraformaldehyde for measuring muscle fiber diameter by hematoxylin-eosin (HE) staining, intramuscular fat content by oil red O staining, and muscle TG, T-CHO content by Tissue TG and T-CHO ELISA kits (Sinobestbio, Shanghai, China) following manufacturer's introductions. A total of 40 one-day-old male Lushi and AA chickens (20 individuals for each breed) were used to measure the birth weight. Then, based on the high uniformity of birth weight, three one-day-old chickens each breed were selected to collect breast muscle samples and equivalently mixed into one composite sample, respectively, for Hi-C analysis.

HE Staining
The breast muscle fiber diameter was measured by hematoxylineosin staining as described by Carvalho (de Carvalho and Taboga, 1996). The paraffin sections of breast muscle were immobilized with 10% paraformaldehyde (Servicebio, Wuhan, China) for 5 min, and were washed with distilled water, then stained with hematoxylin for 15 min. Following 1% hydrochloric acid-ethanol differentiation for several seconds, the paraffin sections were subsequently stained with 0.5% eosin for 3 min and then went through gradient alcohol dehydration for 2 min and vitrification by dimethylbenzene for 2 min. The cell nucleus was stained blue, and the cytoplasm was stained pink. All reactions were performed in triplicate.

Oil Red O Staining
The breast muscle samples (n 3 each group) were sectioned with a freezing microtome. The frozen sections of breast muscle were immobilized with 10% paraformaldehyde (Servicebio, Wuhan, China) for 30 min, washed three times with 1× phosphate buffer saline (PBS) (Gibco, Gaithersburg, MD, United States), and then incubated with oil red O solution (Servicebio, Wuhan, China) for 15 min. Finally, the cell nuclei of breast muscle tissue were restained with hematoxylin (Servicebio, Wuhan, China) for 5 min. The lipid droplet was stained red, and the cell nucleus was stained blue.

Preparation and Sequencing of Hi-C Libraries
The collected fresh muscle samples were cut into pieces with 2-3 mm diameter on the ice and washed in PBS (Gibco, United States). The samples were centrifugated at 4°C, 2,000 × g for 10 min and were resuspended with serum-free Dulbecco's Modified Eagle's medium (DMEM) (Gibco, Gaithersburg, MD, United States). The resuspended cells were DNA-protein crosslinked by a final 2% concentration of formaldehyde, incubated at room temperature for 10 min with soft shock at an interval of 2 min. The cross-linked samples were then supplemented with a final concertation of 0.2 M glycine and incubated at room temperature for 5 min, followed by 15 min on the ice to completely terminate cross-linking, and collected by centrifugation at 4°C, 2,000 × g for 10 min. Total DNA, gelatinin protein residues, and chromatin integrity were used to evaluate the crosslinking effect using the proteinase K DNA extraction method, in concert with agarose gel electrophoresis. The qualified samples were stored at −80°C for further use.
Hi-C experiment was performed as previously reported with minor modifications (Lieberman-Aiden et al., 2009). Briefly, chromatin digestion and isolation were performed by restriction enzyme MboI at 37°C, and the resulting sticky ends were filled in with nucleotides and introduced biotinylated nucleotides at ligation junctions. After that DNA fragment ends were ligated with T4 DNA ligase, purified with a ranging size from 300 to 500 bp, and sheared. Biotinylated junctions were isolated with streptavidin beads, then sequencing connectors were added to form joint products to construct the Hi-C library with an Illumina TruSeq DNA Sample Prep Kit according to the manufacturer's instructions. Hi-C libraries were sequenced on an Illumina Novaseq 6000 platform configured for 150 bp paired-end reads. The Hi-C experiment was performed by Annoroad Gene Technology Co., Ltd. (Beijing, China).

Hi-C Data Analysis
Clean Hi-C data were generated by filtering out low-quality reads (>15% of bases with Q scores ≤ 19%), reads containing over 5% poly-N, and adapters from Hi-C raw data using Trimmomatic (Bolger et al., 2014). Clean reads were mapped to Gallus gallus 5.0 reference genome with a two-step approach embedded in HiC-Pro V2.7.8 software. The unique mapped paired-end reads were generated by filtering out the unmapped paired-end reads, paired-end reads with singleton and multi mapped paired-end reads from clean reads. Then the reads with single side mapping, self-circle ligation, dangling ends, dumped pairs and duplicated valid pairs from PCR artifacts were discarded. Only valid read pairs with unique mapped paired-ends involving two kinds of restriction fragments were used to build the contact maps. The valid paired-end reads were divided into intrachromosomal interactions (cis) reads and inter-chromosomal interactions (trans) reads, where cis rate and trans rate respectively referred to the proportion of cis and trans interaction pairedend reads in the valid paired-end reads. The raw contact maps were normalized with a sparse-based implementation of the iterative correction method by HiC-Pro and visualized by HiTC (Servant et al., 2012). Normalized interchromosomal interaction matrix was obtained from observed/expected number of contacts between all pairs of whole chromosomes to evaluate the spatial proximity of chromosomes. We performed an overall Hi-C subtraction of the Z score matrices (AA minus LS) with genomic distance using a modified LOWESS method, by which the weighted-average and weighted-standard deviation for every genomic distance were calculated to consequently normalize for genomic distance signal bias (Sanyal et al., 2012).

A/B Compartment Analysis
The A/B compartments were estimated by eigenvector analysis of the genome contact matrix after normalization using the observed/expected method (Lieberman-Aiden et al., 2009). Pearson correlation matrix responding to the correlation between each of two loci on the chromosome and the whole genome interaction were employed to dramatically sharpen the plaid pattern. In the normalized interaction matrix heat map, blue indicates a lower observed interaction frequency than expected interaction frequency, on the contrary to red. The correlation coefficient ranged from −1 to 1, and the color transited from blue to red. According to principal component analysis (PCA) analysis with dimensionality reduction of the genome contact correlation coefficient matrix, the chromatin region could be distinctly divided into two parts according to the first eigenvector, which was called the A/B compartment. The positive values and negative values in the first eigenvector respectively, determined the A and B compartments.

TAD Calling and TAD Boundary Analysis
To determine TADs and their boundaries, we proceeded with the previously reported insulation score method (Crane et al., 2015). Briefly, an "insulation score" that reflects the aggregate of interactions occurring across each interval was assigned to genomic intervals along the chromosome to quantify TADs. Of which, the sum of interaction strength of the local chromosome regions was calculated, after smooth transformation, these valleys/minima were interpreted as TAD boundaries or areas of high local insulation, wherein that the Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 782268 minima of insulation profile in total local interaction located at the TAD boundary. It involves two steps, first, the interaction strength value of each bin was calculated by summing its interaction strength in the 10 kb binned Hi-C data (the number of unique mapped read pairs to these regions), subsequently, the interaction strength of each bin served as a value to draw its fluctuation curve on the chromosome to define the boundary.

Prediction and Functional Analysis of miRNAs Target Genes
The miRanda (http://www.microrna.org/microrna/home.do) and TargetScan (http://www.targetscan.org/vert_72/) were employed to predict the potential target genes of miRNA. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of genes was conducted using DAVID (http://david.abcc.ncifcrf.gov/) with a cut-off criterion of p value ≤ 0.05.

RNA Extraction and Complementary DNA (cDNA) Synthesis
Total RNA was extracted from the breast muscle tissues of E10, E14, E18, 1D, 3W, and 5W AA and LS chickens (n 6 each group) using TRIzol ® reagents (Invitrogen, Carlsbad, CA) according to the manufacturer's protocol. The RNA integrity was detected using 1% TAE agarose gel electrophoresis and Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, United States). The concentration, as well as purity assessed by OD260nm/280 nm and OD260nm/230nm, of RNA samples, were measured using NanoPhotometer ® spectrophotometer (IMPLEN, CA, United States). The RNA samples with the ratios of 28S/18S band brightness more than 1.5, as well as 1.8-2.0 of OD260 nm/280 nm and 2.0-2.3 of OD260 nm/230 nm were reversely transcribed into cDNA using PrimerScriptTM RT reagent kit (Takara, Kyoto, Japan) following the manufacturer's instructions.

Statistical Analysis
All data were presented as mean ± standard deviation (SD), and the statistical significance between the two experimental groups was evaluated by t-test for comparisons using SPSS 23.0 (IBM, Chicago, IL, United States). The *p value < 0.05 was considered statistically significant, and **p value < 0.01 was considered extremely significant. Graphics were produced via GraphPad Prism 8 (GraphPad Software, San Diego, CA, United States).

RESULTS
The Differences of Live Weight, Body Size Indexes, Serum Biochemical Indexes, Breast Muscle Fiber Diameter, and IMF Content Between AA Broiler and Lushi Chicken Compared with Lushi chicken, AA broiler is characterized by extremely significant higher live weight (p 3.35E-22) and body size indexes at 5 weeks old including shank length (p 3.68E-19), shank circumference (p 1.47E-14), sternal length (p 2.02E-22), body length (p 8.87E-27), chest depth (p 3.23E-23), chest width (p 9.19E-23) and chest angle (p 3.31E-20) ( Table 1). There was no significant difference in the serum levels of ALT (p 0.581), AST (p 0.367), GLU (p 0.083), LDH (p 0.823), and TG (p 0.112) between 5-week-old AA broiler and Lushi chicken, although serum TG exhibited an increased trend in LS chicken. Of note, significantly higher serum T-CHO (p 0.005), HDL-c (p 0.005), LDL-c (p 0.001), and VLDL-c (p 0.047) were found in LS chicken versus AA broiler ( Table 1). The remarkably increased IMF content of pectorals ( Figure 1A), but remarkably decreased breast muscle fiber diameter ( Figure 1B) were observed in Lushi chicken in comparison with AA broiler at 5 weeks old. Simultaneously, the conspicuously superior IMF content was uncovered in one-day-old LS chicken compared to the AA broiler ( Figure 1C), but the latter exhibited a prominent breast muscle fiber diameter in pectorals ( Figure 1D). Moreover, significantly increased TG level was observed in the muscle of one-day-old LS chicken compared to AA (p 0.0010) ( Figure 1E), and there was no significant difference of muscle T-CHO between them (p 0.1048) ( Figure 1F).

Characterization of Hi-C Data From AA and LS Breast Muscles
As shown in Table 2

Comprehensive Maps of Chromatin Interaction in the Breast Muscles of AA and Lushi Chickens
Whole genome-wide intra-(cis) and inter-(trans) interaction data of AA and LS chicken breast muscles were visualized as chromosome versus chromosome heat maps, where darker colors represent more frequent interaction events (Figures 2A,B). Of which, 89.98% of AA interaction data and 83.4% of LS interaction data were found as cis interaction, 10.02 and 16.6% as trans interaction ( Figure 2C). Then, to evaluate whether the chromosomes clustering was altered between AA and LS chicken breast muscles, we compared the chr1-20 and chrZ genome-wide interaction differences and draw LS minus AA genome-wide heatmap of significant differential interactions. It was demonstrated that compared with AA, LS was equipped with significant increases in genome-wide interaction ( Figure 2D). In order to intuitively show the trans-interaction, 1,000 bin pairs with the strongest interaction among chromosomes were used to draw circle plots. The circle plots showed that there was frequent interaction between the small, gene-rich, and high gene density chromosomes in AA and LS chicken breast muscles. Compared to LS, AA was equipped with the relatively repressed interaction frequency between large chromosomes and small chromosomes, in parallel with more frequent interaction between the large chromosomes ( Figures 2E,F). Moreover, we standardized the interaction between each chromosome using the observed/ expected number of contacts between all pairs of whole chromosomes to get a standardized chromosome interaction matrix in AA and LS chicken breast muscle, where the higher interaction between two chromosomes represents closer physical proximity ( Figures 2G,H). As such, AA broiler breast muscle was characterized by more frequent interaction between small, gene-rich chromosomes, especially Chr25-33, in parallel with LS showing a similar trend except for the weaker interaction between Chr16 and Chr17-33. Considering that chromosome size made a difference in chromosomal interaction, we sorted all chromosomes subjected to gene density and conducted hierarchical clustering to exhibit the correlation more intuitively between chromosome size and the intensity of chromosomal interaction according to chromosomal interaction values. It was revealed that Chr30 with a size of 0.025 Mb possessed the highest gene density, following Chr32, Chr31, Chr16, Chr25, Chr33, Chr28, Chr27, Chr26, Chr23, Chr21 raging in size from 0.04 to 7 Mb ( Figure 2I). Interestingly, hierarchical clustering analysis of whole chromosome positioning indicated that the chromosome trans-interaction was heavily concentrated on the small and gene-rich chromosomes in AA and LS chicken breast muscle ( Figures 2J,K).

Identification and Switching Characterization of A/B Compartments in the Breast Muscles
To verify whether there are any differences in the compartmentalization between the AA and LS genomes from breast muscles, we compared the A/B compartments throughout the whole genome at 100 kb resolution. There was no significant difference in the number of B type compartments between the two breeds, except more A type compartments in AA versus LS ( Table 3). AA and LS chicken breast muscles shared 47.69% stable A to A compartments and 48.73% stable B to B compartments, where 2.02% of all compartments constituted a switching in genomic compartmentalization from A-type in AA to B-type in LS, referring to A to B switching, and 1.56% vice versa ( Figure 3A; Supplementary Table S1). To further explore whether small and gene-rich chromosomes dominatingly contributed to the A/B compartment switching, the compartment distribution of Chr16-33 and the rest of the genomes were analyzed. A distinct increase of stable A compartments and compartments switching from A-type in AA to B-type in LS on these chromosomes were found ( Figure 3B). The small, gene-rich chromosomes possessed relatively abundant A compartment and poor B compartment  Figure S1). We emphasized the A/B compartment switching on Chr21-33 and observed that there was no compartment switching on Chr21, on which both AA and LS merely constituted 1.45% B to A compartment switching, together with 52.17% stable A compartments and 46.38% stable B compartments ( Figure 3C; Supplementary Figure S2A Figure 3F; Supplementary Figure S2D). In terms of Chr28, both AA and LS merely constituted stable A compartments and 44% stable B compartments ( Figure 3G; Supplementary Figure  S2E). Only A to B compartment switching was respectively, detected on Chr24 (3.17%), Chr25 (2.86%), Chr27 (1.75%), and Chr33 (11.76%) in the breast muscles of the two kinds of chickens ( Figures 3H-K Figure 4A; Supplementary Table S2).
To better understand the correlation of genomic compartmentalization with regulating muscle development and IMF deposition between AA and LS, we performed the functional annotation and signaling pathway enrichment analysis of genes located in the A/B switching region in their breast muscles. A total of 25 clusters based on the GO functional annotation of mRNAs in the A to B compartment switching region were validated. There were nine significant terms including negative regulation of platelet aggregation (p 0.0030), cell proliferation (p 0.0132), blood coagulation (p 0.0253), positive regulation of neuron projection development (p 0.0253), feeding behavior (p 0.0395), negative chemotaxis (p 0.0436), an integral component of the plasma membrane (p 0.0123), thrombin receptor activity (p 0.0019), and transcriptional activator activity, RNA polymerase II distal enhancer sequence-specific binding (p 0.0224) ( Figure 4B; Supplementary Table S3). Moreover, 19 clusters based on the GO functional annotation of mRNAs in B to A compartment switching region were validated, including five significant terms: cellular calcium ion homeostasis (p 0.0334), dystroglycan binding (p 0.0016), peptidyldipeptidase activity (p 0.0382), SNAP receptor activity (p 0.0469) as well as transcription, DNA-templated (p 0.0351) ( Figure 4C; Supplementary Table S3). Moreover, the top 20 GO terms from mRNAs in stable A and B compartments according to p value were respectively, selected and showed that the majority of terms focused on cellular components, following molecular function and biological process least in both stable A and B compartments (Supplementary  Figures  S3A,B  and  Supplementary Table S3).
To identify the crucial signaling pathways that participated in muscle development and IMF deposition between AA and LS, we mapped these genes respectively, from stable and switched compartments to KEGG orthologs. The genes on A to B compartment switching region were significantly enriched in the neuroactive ligand-receptor interaction pathway (p 0.0129). Genes in the B to A compartment switching region were enriched in the calcium signaling pathway (p 0.0632) and TGF-β signaling pathway (p 0.0761). Concerning mRNAs on stable A compartment, they were significantly enriched in the biosynthesis of antibiotics (p 2.15E-05), pyruvate metabolism (p 0.0015), insulin signaling pathway (p 0.0052), carbon metabolism (p 0.0055), protein processing in the endoplasmic reticulum (p 0.0095), spliceosome (p 0.0228), cell cycle (p 0.0236), mRNA surveillance pathway (p 0.0344), regulation of autophagy (p 0.0384), nucleotide excision repair (p 0.0410), purine metabolism (p 0.0459) and two amino acid-related metabolism pathways including glycine, serine and threonine metabolism (p 0.0097), and biosynthesis of amino acids (p 0.0254). The mRNAs on the stable B compartment were significantly associated with neuroactive ligand-receptor interaction (p 3.96E-15), cell adhesion molecules (CAMs) Given that miRNA can act as a post-transcriptional regulator to participate in various physiological processes via complementarily binding to the 3′-untranslated region (3′UTR) of its target, to comprehensively investigate the potential effects of A/B compartment switching on muscle development and IMF deposition, the targets of miRNA from A to B compartment switching region and B to A compartment switching region were analyzed, respectively. A total of 4,722 genes were obtained from two intersecting software programs, revealing that miRanda and TargetScan were potential targets of miRNAs located on the A to B Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 782268 compartment switching region. They were mainly clustered in regulatory and metabolic processes, otherwise directly related terms to lipid metabolisms such as acyl-CoA dehydrogenase activity, coenzyme A biosynthetic process lipid homeostasis, and fatty-acyl-CoA binding ( Figure 4E; Supplementary  Table S5). Moreover, 1,479 potential targets of miRNAs located on the A to B compartment switching region were predicted and clustered in cellular activity, mRNA processing, and organelle-related terms ( Figure 4F; Supplementary Table S5). The signaling pathway enrichment analysis indicated that the potential targets of miRNAs were significantly associated with fatty acid metabolism (p 0.0029), protein processing in the endoplasmic reticulum (p 0.0092), cell cycle (p 0.0109), biosynthesis of unsaturated fatty acids (p 0.0162), fatty acid elongation (p 0.0162), glutathione metabolism (p 0.0173), valine, leucine, and isoleucine degradation (p 0.0352), pyrimidine metabolism (p 0.0497) on the A to B compartment switching region, as well as endocytosis (p 2.74E-05), apoptosis (p 3.63E-05), protein processing in the endoplasmic reticulum (p 0.0008), cell cycle (p 0.0048), FoxO signaling pathway (p 0.0115), influenza A (p 0.0406) on the B to A compartment switching region ( Figure 4G; Supplementary  Table S5).

Topologically Associated Domains Calling and Boundary Sliding in Breast Muscle Between the Two Breeds
The TAD boundaries were responsible for maintaining TAD independence and structural integrity, which are necessary for gene regulation. To investigate whether TAD boundaries were altered by large-scale chromosomal interactions and genomic compartment switching between AA and LS chicken breast muscle genomes, the TAD boundaries were identified by calculating the insulation plot of genome-wide interaction maps at 40 kb resolution and found that 357 and 356 TAD boundaries in AA and LS genomes, respectively. Among those, there were 178 overlapping TAD boundaries in both AA and LS, 179 AA-specific TAD boundaries, and 178 LS-specific TAD boundaries ( Figure 5A). Figure 5B, AA and LS did not show a significant change in TAD score distribution on whole chromosomes. In terms of the average size of TADs, there was no significant change between AA and LS ( Figure 5C), wherein TADs on Chr10 were charactered by the largest average size in length in both AA and LS genomes ( Figures 5D,E). The TAD distributions on all chromosomes of AA and LS were shown in Figures 5F,G, respectively, indicating Chr1 possessed a maximum number of TAD and Chr17 possessed a minimum in these chromosomes with TADs. It has been well acknowledged that the dynamic changes of TAD boundaries give rise to the reforming and/or disappearing of TAD, resulting in the alternation of regulatory activity and following abnormal gene expression and subsequent disease phenotypes. There were 29 reforming events (Chr1,2,3,4,6,7,8,14,20,Z) and 32 disappearing events (Chr1,2,3,4,5,6,13,15,Z) of TAD boundaries existing in AA versus LS chicken breast muscle. A total of 50 and 67 genes, were located on AA and LS TAD boundaries, respectively (Supplementary Table S6). And the 50 genes were clustered into two GO terms including defense response to the bacterium (p 0.00402) and extracellular region (p 0.00509). The 67 genes were clustered into 11 terms, including defense response to Gram-negative bacterium (p 1.2575E-05), defense response to Gram-positive bacterium (p 1.41986E-05), innate immune response (p 0.00086), membrane disruption in other organisms (p 0.00378), negative regulation of gene expression (p 0.00672), mitochondrion (p 0.00891), cytolysis (p 0.02061), positive regulation of gene expression (p 0.02377), lipopolysaccharide binding (p 0.02779), developmental growth (p 0.03351) as well as cell death (p 0.03717) ( Table 4). Of which, insulin like growth factor 2 mRNA binding protein 3 (IGF2BP3), a histone modification factor as m6A readers involving m6A RNA methylation were found in AA-specific TAD boundary, implying that there was probably a connection with the formation of TAD ( Figure 6A). Another lipid-related gene, 3-hydroxy-3-methylglutaryl-CoA reductase (HMGCR) was located in the AA-TAD boundary, whereas in LS-TAD interior regions ( Figure 6B). To determine the difference in the mRNA abundance of the candidates, the relative expression levels of IGF2BP3 and HMGCR genes in each developmental stage of chicken breast muscle were quantified by qRT-PCR. The expression of the IGF2BP3 gene was gradually decreased as breast muscle development was processed in both AA and LS chickens. However, a significantly increased IGF2BP3 mRNA expression was observed in the breast muscle of the AA broiler compared to LS chicken during all breast muscle development stages (p < 0.05) ( Figure 7A). Similarly, the mRNA abundance of the HMGCR gene exhibited a gradually decreased trend during breast muscle development stages in LS chicken and AA broiler, while HMGCR mRNA expression was overall up-regulated in breast muscles of AA broiler than LS chicken, and was especially significantly higher at E10 (p 0.049), E14 (p 0.014), E18 (p 0.020), 1D (p 0.020), and 5W (p 0.001) ( Figure 7B).

DISCUSSION
Chicken has been widely accepted as a crucial research model organism for implications on phylogenetics and embryology, along with agriculture and biology, such as adipogenesis, development, immune, medicine, and cancer, owing to its convenience for operation and the extensive conservation of synteny between chicken and mammals (Burt, 2007;Beacon and Davie, 2021). Both muscle development and IMF deposition are complex and precisely orchestrated processes mediated by a network of regulatory factors directly and indirectly related to growth development and lipid metabolism. Increasing investigations are committed to elucidating these regulatory mechanisms at genomic, transcriptomic, and proteomic scales, where the genome has been deemed a linear molecular model. However, this is insufficient when it comes to revealing the connection between discrete regulatory elements, structural variations, and gene function, and consequently, it is imperative to apply threedimensional genomics based on the spatial conformation of chromatin to explain the regulatory mechanism of gene expression. Therefore, in this study, we first identified a comparative analysis of genome-wide interactions between AA and LS chicken breast muscles upon hatching, when IMF has been proved to accumulate earliest, using Hi-C to illuminate the regulatory mechanism underlying chicken muscle development and IMF deposition at the three-dimensional chromosome scale. Different from mammals whose lipid metabolism occurs in adipose tissue, it is the liver that principally governs lipid metabolism in chicken, where over 90% de novo fatty acids are synthesized. Endogenous lipids synthesized by the liver and exogenous sources absorbed from the diet undergo catalytic reaction into acyl-CoA and subsequently metabolization to form TG, phospholipids, and total cholesterol, which are incorporated into very low-density lipoprotein (VLDL) by apolipoprotein and then transported in the blood circulation (Leveille et al., 1968;O'Hea and Leveille, 1969). After that, VLDL is hydrolyzed into fatty acids by lipoprotein lipase (LPL) and then uptaken by fat cells in muscles, and reesterified to form TG, which is stored to enhance the IMF deposition (Qiu et al., 2017). Therefore, serum TG, T-CHO, and VLDL content mostly reflect lipid metabolism activity in chicken, and serum VLDL has been widely regarded to have a strong correlation with fat deposition (Griffin and Whitehead, 1982;Hermier, 1997). Our results showed the increased serum TG, T-CHO, and VLDL levels, together with enhanced accumulation of lipid droplets in breast muscle of LS chicken compared with that in breast muscle of AA broiler. It implied that more IMF deposition in LS chicken breast muscle might benefit from dramatically increased lipid metabolism in LS compared with AA, consistent with those signaling pathways related to lipid metabolism could contribute to the high capacity of IMF accumulation (Cui et al., 2012;Liu R. et al., 2017). Interestingly, LS chicken has markedly superior IMF content to AA broilers, while there was no significant difference of muscle T-CHO between them, suggesting a lower intramuscular cholesterol proportion of fat in breast muscle of LS chicken than that of AA broiler, indicating that the meat has lower cholesterol. Moreover, the accretion rate of IMF depends upon the muscle growth rate, that is, increasing muscularity could dilute the final fat content of muscle (Hocquette et al., 2010). Here, our study showed a larger breast muscle fiber diameter in AA than that in LS, to some extent, those might also explain why more IMF deposition existed in LS chicken breast muscle. As expected, the normalized chicken Hi-C maps of AA and LS chicken breast muscle displayed a typical plaid-pattern and strong signals along diagonals, rooted in striking nonrandom patterns in its genomic sequence composition (Bickmore and van Steensel, 2013). The interactions occurred most frequently in cis and decay with genomic distance, in agreement with a previous report revealing that the intensity of chromatin interaction decreased as the linear distance of the genome increased, and is higher within the same chromosome than between different chromosomes (Lieberman-Aiden et al., 2009). In mammalian cell nuclei, the distribution of chromosomes is not random but closely concerned with gene density, chromosome size, GC content, and transcription activity, etc, wherein the gene-rich chromosomes concentrate at the center of the nucleus and the gene-poor chromosomes are located towards the nuclear periphery (Sun et al., 2000;Boyle et al., 2001;Parada et al., 2003;Misteli, 2007;Lieberman-Aiden et al., 2009). Our studies have revealed that, in both LS and AA breast muscle, the small gene-rich chromosomes displayed increased interchromosomal interaction frequency between pairs of chromosomes compared with large gene-poor chromosomes, corresponding to these chromosomes being equipped with a widely acknowledged strong physical proximity and frequent colocalization in the center of the nucleus. The formation and disappearance of chromatin interactions have been widely acknowledged to play a critical role in cell fate determination and differentiation, cooperating with gene-specific expression regulation (Bartkuhn and Renkawitz, 2008). Of note, interchromosomal interaction frequency showed that the large, gene-poor chromosomes in LS chicken breast muscle preferentially interacted with each other compared to AA broiler, speculating that it might contribute to explaining  the difference of muscle development and IMF deposition, which needs to be further studied.
It has been acknowledged that compartment A has a more open spatial structure, higher gene density, stronger gene expression activity, and higher GC base content; on the contrary, chromatin compartment B contains centromeres and has a more closed spatial structure with lower internal gene density and lower transcriptional activity. Within the same type of compartment, chromatin interaction frequency is higher, and at the same linear distance, the interaction frequency between B type compartments is superior to A type compartments. In our study, the B compartment of either LS or AA chicken breast muscle covered approximately 60% of the chicken genome, in accordance with the point in rice (Liu C. et al., 2017). There was no significant difference in the number of B type compartments between the two breeds, more A type compartments in AA being likely to cause the dynamic regulation of gene expression related to muscle growth and development and/or lipid metabolism, resulting in a difference in muscle development and IMF deposition between AA and LS chickens. During muscle development, the programmed, controllable and orderly expression of a large number of genes leads to the specific expression of muscle cell genes, forming a complex regulatory network and signaling pathways. Simultaneously, IMF deposition is also mediated by multiple pathways and influenced by a number of genes, especially in chickens, because, unlike mammals, the liver is principally responsible for lipid metabolism in avian. According to our gene annotation analysis of mRNAs located in stable A and B compartment regions, these genes were involved in the development, enzymatic activity, transport, immune response, and other processes. The chromatin organization and stability are considered to be regulated by post-translationally histone modifications, such as methylation, acetylation, phosphorylation, and ubiquitination (Elgin and Workman, 2002;Shah, 2015). Moreover, there is a mutual dependency between chromatin organization and patterns of epigenetic marks. For example, DNA methylation is confined to the A compartment in cardiac myocytes (Nothjunge et al., 2017), unmethylated CpGs are enriched in the A compartment, and methylation levels are decreased to a greater extent in the A compartment than in the B compartment in mouse embryos (Ke et al., 2017). Moreover, both global demethylation and remethylation correlate with chromatin compartments in early embryonic development in mammals (Zhang et al., 2018), etc. Interestingly, within A compartments, there are some processes involving histone modifications, for instance, histone acetyltransferase activity, ubiquitin protein ligase activity, ubiquitin protein ligase binding, ubiquitin conjugating enzyme activity, protein dephosphorylation, ubiquitin-dependent protein catabolic process, protein autoubiquitination, ubiquinone biosynthetic process, histone H4 acetylation, suggesting that more histone modifications take place in the A compartment in chicken breast muscle. Moreover, the genes located within A compartments participate in various pathways related to energy metabolism, autophagy, fatty acid metabolism, insulin signaling pathway, and Notch signaling pathway, which have been demonstrated to be involved in muscle development and/or lipid accumulation. For example, Lactate Dehydrogenase A (LDHA), acyl-CoA synthetases (ACSs) gene family, long chain fatty acid elongases (ELOVLs) gene family, stearoyl-CoA desaturase (SCD) and histone deacetylases (HDACs) gene family, hydroxyacyl-CoA dehydrogenase (trifunctional protein), alpha subunit (HADHA), HMGCR, and mitogenactivated protein kinases (MAPKs) gene family are universal genes responsible for the development of skeletal muscle and lipid metabolism (Mashek, 2013;Li et al., 2016;Kim et al., 2017). From the late embryonic stage to hatching, a great number of genes also undergoe transcription and translation to form proteins to meet the following demands of energy at hatching. This could be the reason that the genes concerned with energy metabolism are located in the A compartment in the chicken breast muscle at hatching. The alterations of compartment types (from A to B, or from B to A) mean that changes of chromatin state (from an open active state to a closed inactive state, vice versa) directly affect the inhibition or activation of gene expression. Here, the KEGG pathways showed that the genes within A to B compartment switching regions were clustered to Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 782268 neuroactive ligand-receptor interaction, which has been proven to exert a role in adipocyte differentiation in human adiposederived stem cells (Guo and Cao, 2019) and bovine muscle development (Cassar-Malek et al., 2007). Furthermore, the genes within the B to A compartment switching regions were significantly involved in the calcium signaling pathway, a highly versatile intracellular signal being capable of regulating many different processes (Berridge, 2012), as well as the TGF-β signaling pathway, which could multifunctionally regulate many cellular processes including cell growth, cell differentiation, apoptosis, morphogenesis, immune regulation, wound healing, inflammation, and cancer in both adult organisms and developing embryos (Morikawa et al., 2016).
These results indicate that A/B compartment switching might yield the dynamic regulation of biological processes, leading to the difference of muscle development and IMF deposition between LS and AA chicken breast muscle. It is of note that the representative genes related to lipid metabolism, such as monoacylglycerol O-acyltransferase 1 (MOGAT1), acyl-CoA synthetase long chain family member 3 (ACSL3), lipin 1 (LPIN1), CCAAT/enhancer binding protein gamma (CEBPG), together with the representative genes related to muscle growth and development, such as insulin like growth factor 2 mRNA binding protein 1 (IGF2BP1), insulin like growth factor 1 (IGF1), myosin, and heavy chain (MYH) gene family were found within the A to B compartment regions. Regarding B to A compartment switching regions, representative genes including Wnt family member 16 (WNT16), SET and MYND domain containing 2 (SMYD2), BTG anti-proliferation factor 2 (BTG2), and PR/SET domain 6 (PRDM6) have been confirmed to exert vital functions in both muscle development and lipid metabolism in mammals. The A/B compartment switching might generate the changes of expression levels of these genes, thereby forming the distinct ability of muscle development and IMF deposition between LS and AA chickens. However, further research is needed to identify which are the dominant genes. An increasing number of investigations have indicated that TAD is not only a structural unit of chromatin but also a functional unit involved in many biological processes, such as genome evolution, DNA replication, and gene regulation (Galupa and Heard, 2017;Szabo et al., 2019). The stability of the TAD structure plays a crucial role in maintaining the systematical and sequential physiological activities of the organism. The reconstruction and/or disappearance of TAD boundaries would directly lead to abnormal gene expression, consequently causing diseases such as human finger deformities (Lupiáñez et al., 2015) and fetal severe limb anomalies (Trimouille et al., 2019). We found no significant difference in TAD distribution and size occurred in LS and AA chicken breast muscles. Almost 50% of all the TAD boundaries between LS and AA chicken breast muscles were overlapped, indicating that the TAD positioning in breast muscles between the two breeds was moderately conserved, and that, their transcriptional expression regulation varied. This might partly be limited by sequencing technology and the quality of the chicken reference genome. Similar to a previous report, it is hard to identify TADs within several of the smallest chicken chromosomes, which are poorly represented in the current assembly of the chicken genome (Schmid et al., 2015). While further improvement of the chicken genome assembly emerges, the TAD analysis of these tiny chromosomes in the interphase nucleus could be supplemented. It has been verified that the TAD boundary was enriched with histone modification, transcriptional start sites, RNA polymerase II, and transcriptional factors, and could be associated with the different gene distribution and expression pattern from TAD-interior regions (Lieberman-Aiden et al., 2009;Mourad and Cuvier, 2015;Hong and Kim, 2017). Consistent with this opinion, our study also showed that TAD boundaries comprised histone modification factors and transcriptional factors, indicating a certain connection between transcription and local chromatin topology in AA and LS chicken breast muscles.
It is noteworthy that the IGF2BP3 gene, a family of RNAbinding proteins that also acts as an RNA N6-methyladenosine reader, is located on the AA unique TAD boundaries, which were related to mRNA stability and translation. It could regulate the production, localization, and expression of downstream genes, especially IGF2, acting as a master switch governing the initiation of skeletal muscle development (Nielsen et al., 1999;Chao and D'Amore, 2008;Huang et al., 2018). A recent study demonstrated that knockdown of the IGF2BP3 gene significantly delayed differentiation and induced proliferation of C2C12 myoblasts, and revealed that IGF2BP3 might be considered as a candidate gene in pig breeding for meat production traits (Yang et al., 2021). Additionally, IGF2BP3 has been identified as a potential candidate gene related to IMF deposition in chicken, and the mRNA expression of IGF2BP3 gene in the insulin pathway was decreased in skeletal muscle of dwarf chicken, whose IMF content was significantly increased, compared with normal chicken (Lin et al., 2012;Ye et al., 2014), which is consistent with our present study that there is lower IGF2BP3 expression in LS chicken than fast-growing AA broiler. This evidence further implies that IGF2BP3 might play indispensable roles in breast muscle development and intramuscular fat deposition in chicken. HMGCR could drive the catalyzed conversion of HMG-CoA into mevalonic acid, a rate-limiting step in cholesterol synthesis, thus playing a critical role in cholesterol homeostasis (Sharpe and Brown, 2013). It has been reported that the polymorphisms of the HMGCR gene were associated with several lipid deposition-and cholesterolrelated traits, and also with IMF content of gluteus medius muscle in pigs (Canovas et al., 2010;Chalupová et al., 2012). In bovine intramuscular adipocytes, the alteration of the HMGCR gene could affect the expression abundance of AMP-activated protein kinase (AMPK) and sirtuin type 1 (SIRT1), which exert important roles in the regulation of energy metabolism (Lin et al., 2020). The polymorphisms in the 3′UTR of the HMGCR gene exhibited significant association with the total cholesterol content in chicken muscle. This evidence shows a strong correlation between HMGCR gene expression with IMF deposition in livestock and poultry (Cui et al., 2010). In our present study, the decreased mRNA expression of the HMGCR gene was detected in the breast muscle of LS chicken compared to the AA broiler, in accordance with the lower T-CHO proportion of LS chicken, Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 782268 leading to the speculation that LS could be managed by a slower rate of intramuscular cholesterol synthesis. It was also evidenced that anti-HMGCR could inhibit muscle cell fusion and trigger atrophy on fully differentiated myotubes in vitro primary human myoblasts/myotubes isolated from human muscle biopsies of nonmyopathic patients, suggesting that HMGCR could facilitate muscle development in humans (Arouche-Delaperche et al., 2015). Additionally, RNAi-HMGCR specifically in the corpus allatum yields dwarf flies, suggesting HMGCR could control development in Drosophila (Belgacem and Martin, 2007). In chicken, polymorphism of the HMGCR gene was reported to correlate with leg muscle weight and leg muscle fiber diameter, revealing that the HMGCR gene might also exert a crucial role in muscle development (Wei et al., 2012). The lower levels of IGF2BP3 and HMGCR mRNA expression during breast muscle development in LS chicken might contribute to the suppressive effect on breast muscle development and the stimulative effect on IMF deposition. Given that classical TAD-boundary regions were linked to the transcriptional control of mammalian genomes (Dixon et al., 2012), it has been proven that, compared to TADinterior regions, the genes located in TAD-boundary regions were prone to being highly expressed (Dong et al., 2018). This might explain the higher expression profiles of HMGCR genes in one-day-old AA chicken breast muscles, which are located in the AA TAD-boundary regions, but in LS TAD-interior regions, as well as higher expression profiles of IGF2BP3 gene located on the AA unique TAD boundaries. Collectively, this implies that the TAD boundary disappearing or sliding in LS chicken verse AA chicken could contribute to lower transcriptional expression of IGF2BP3 and HMGCR genes in LS chicken breast muscles, and was thus responsible for slower muscle development as well as higher IMF content and lower intramuscular cholesterol proportion in LS chicken, thereby suppressing the breast muscle development but yielding excellent meat quality in LS chicken. However, the specific biological functions and regulatory mechanism at chromatin level underlying chicken breast muscle development and IMF deposition require further research.

CONCLUSION
In conclusion, the chromatin structure of the chicken breast muscle from LS and AA broiler chickens was identified and compared using large-scale chromosomal cis-and transinteractions to examine genomic compartmentalization (A/ B compartments) and TAD formation. Two genes IGF2BP3 and HMGCR regulated by TAD boundary sliding were identified as potential biomarkers for chicken breast muscle development and IMF deposition. To our best knowledge, our findings are the first to illuminate the regulatory mechanism of muscle development and IMF deposition at a threedimensional chromosomal level in chickens, providing new insight into the functional role of chromatin organization in the formation of avian economic traits.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/, PRJNA757147.

ETHICS STATEMENT
The animal study was reviewed and approved by the Institutional Animal Care and Use Committee of Henan Agricultural University (Permit Number: 11-0085).

AUTHOR CONTRIBUTIONS
WT performed the experiments, analyzed data, and prepared the manuscript. ZW and DW collected the samples and participated in immunohistochemistry study. YZ and JD performed the RNA extraction and qRT-PCR. RJ, RH and ZL participated in data analyses and provided comments on manuscript revision. XK participated in the design of the study and critical discussion on the results. HL participated in the design of the study and helped to revise the manuscript. XL conceived the study and provided overall supervision. All authors read and approved the final version of the manuscript.