Abstract
The spatial organization of the genome in the nucleus plays a crucial role in eukaryotic cell functions, yet little is known about chromatin structure variations during late fetal development in mammals. We performed in situ high-throughput chromosome conformation capture (Hi-C) sequencing of DNA from muscle samples of pig fetuses at two late stages of gestation. Comparative analysis of the resulting Hi-C interaction matrices between both groups showed widespread differences of different types. First, we discovered a complex landscape of stable and group-specific Topologically Associating Domains (TADs). Investigating the nuclear partition of the chromatin into transcriptionally active and inactive compartments, we observed a genome-wide fragmentation of these compartments between 90 and 110 days of gestation. Also, we identified and characterized the distribution of differential cis- and trans-pairwise interactions. In particular, trans-interactions at chromosome extremities revealed a mechanism of telomere clustering further confirmed by 3D Fluorescence in situ Hybridization (FISH). Altogether, we report major variations of the three-dimensional genome conformation during muscle development in pig, involving several levels of chromatin remodeling and structural regulation.
1 Introduction
Deciphering the mechanisms that govern gene expression regulation is essential for understanding the fundamental biological changes occurring under different physiological conditions. In this context, genome organization has been proven to be a major player in the regulation of gene expression (; ; ). Understanding the relationship between genome organization and gene expression needs a deep knowledge of chromatin structure and folding, which has been made possible by the development of three-dimensional (3D) techniques like 3D DNA Fluorescence in situ Hybridization (FISH) and Chromosome Conformation Capture assays (; ), including its genome-wide version Hi-C (). By identifying pairs of genomic regions in direct physical contact or in close spatial proximity within the nucleus, hereafter referred as “interactions”, these approaches revealed several features of the genome architecture. For instance, individual chromosomes occupy discrete territories in the interphase nuclei, the so-called chromosome territories (; ; ), which may intermingle at the interface regions allowing trans-chromosomal interactions (; ). Moreover, chromosomes have been found to be organized in two main types of large regions with different features in terms of genome topology, chromatin state and gene expression. These regions of several megabases are the A and B compartments that correspond respectively to open transcriptionally active and close inactive chromatin. While A compartments are associated with euchromatic, gene-rich and DNase I hypersensitive regions, B compartments are considered as transcriptionally inert, heterochromatic, nuclear lamina-associated, gene-poor and DNase I insensitive (; ; ). Although these compartments could be further segmented considering finer epigenetics features () or associated with exceptional euchromatin/heterochromatin organisations (), we will simply refer to the general A/B definition hereinafter. At a smaller scale, genomic regions of about 1 Mb with a high density of cis-interactions, named Topologically Associating Domains (TADs) (; ; ) have been shown to play a role in regulating gene expression during key biological processes like development (; ).
To gain insight into the establishment, the dynamics and the function of these genomic structures, several studies have characterized them in various cell types and compared them within or, sometimes, between species (; ; ). Various comparisons have been performed during early embryo development (Zheng and Xie, 2019), between different cell lines (e.g., embryonic and mesenchymal stem cells ), from distinct differentiation states (e.g., during neural differentiation , or during B cell fate commitment ; ). Such comparisons efficiently revealed strong differences between distinct cell types, requiring few biological replicates (often simple duplicates), but they did not provide information about the heterogeneity and the dynamics of the genome 3D structure for a specific cell type. The development of single-cell Hi-C () made possible to determine whole genome structures in single haploid () or diploid cells (). Recent applications of single-cell Hi-C revealed various degrees of heterogeneity in genome 3D conformation among several cell lines (; ).
Despite all these efforts, little is known about the status and the dynamics of chromosome organization in animal cells from most of the organized tissue types, with notable exceptions like brain and liver for instance (; ; ; ). Regarding skeletal muscle, Hi-C experiments have been performed on cultured cells (; ) and on adult muscle (), but little is known about chromosome organization in this type of differentiated cells during late development. To assess whether significant structural dynamic modifications could also be detected there, we characterized the 3D genome organization of porcine longissimus dorsi muscle cells during late fetal development (days 90 and 110 of gestation) by adapting the in situ Hi-C protocol () to fetal frozen tissues. This period, which covers almost the entire last month of gestation, is known to be crucial for porcine muscle development and maturity, involving major reorganizations of the transcriptomic and proteomic programs (, ). By performing the experimental assays on tissue samples from different fetuses (three replicates per group of the Large White breed), we characterized the genomic structure of pig muscle cells at various levels of organization, providing high-resolution Hi-C interaction maps, TAD and A/B compartment annotations. Comparing samples from 90 vs 110 days of gestation allowed the identification of major topological differences between the two groups, in line with previous results from transcriptome characterization. In addition, these results completed and further expanded previous studies which identified trans interactions involving genes that are key players for fetal muscle growth and development (; ). Overall, this study sheds a new light on the description of dynamic changes of the 3D genome occurring during transcriptional switches in the expression programs of differentiated cells.
2 Materials and Methods
Experimental and computational resources used in this study are listed in Supplementary Table 3, including names of chemical reagents, kits and software versions.
2.1 Animals and Samples
For Hi-C and FISH experiments, longissimus dorsi fetal porcine muscle samples were collected from the European Large White (LW) breed (F1 LW × LW). Specifically, three 90 days gestation male littermates and three 110 days gestation (two male littermates and one female) were used for Hi-C assays. For FISH experiments, muscle samples were collected from different fetuses (one at 90 days gestation and one at 110 days) of those in which Hi-C experiments were performed. All the fetuses used in this study were obtained by caesarean after euthanasia of healthy wild type sows and fetuses. No special breeding conditions (feeding, housing, treatment) were applied.
The experimental design was approved and authorized by the ethical committee (#84) in animal experimentation of the French Ministry of National Education, Higher Education, and Scientific Research (authorization #02015021016014354). The experiment authorization number for the experimental farm GenESI (Genetics, testing and innovative systems experimental unit) is A17661. The procedures performed in this study and the treatment of animals complied with European Union legislation (Directive 2010/63/EU) and French legislation in the Midi-Pyrénées Region of France (Decree 2001–464). All the details about the animals and samples have been registered in the BioSamples public repository (https://www.ebi.ac.uk/biosamples) in agreement with the FAANG best practices guidelines (https://www.faang.org/data-share-principle) and are available using the accession SAMEA7390788.
2.2 3D DNA FISH Experiments
2.2.1 Cells and Probes Preparation
Fetal muscle tissue was obtained from the Longuissimus dorsi muscle of 90- and 110-days of gestation Large White (LW) pig and prepared as described in ; . Stored muscle fibre packets were permeabilised for 5–8 min in cytoskeleton extraction buffer (100 mM NaCl, 300 mM sucrose, 3 mM MgCl2, 10 mM PIPES pH 6.8) containing 0.5% Triton X 100 and then fixed in cold 4% paraformaldehyde for 5 min. After washing in cold PBS, muscle packets were manually dilacerated directly on Superfrost glass slides (CML, Nemours, France) to isolate individual fibres, and air-dried before adding DNA probes for in situ hybridization. Bacterial artificial clones (BACs) containing specific subtelomeric sequences of porcine chromosomes 2, 9, 13 and 15 were chosen as selected by : SSC2p (PigI-370D12), SSC9q (PigI-441D12, PigI-564B6), SSC13q (PigI-39F7) and SSC15q (PigI-899B10). These BACs were isolated from a porcine BAC library (CRB-Anim, INRA, 2018. Biological Resource Centres for domestic animals of AgroBRC, doi: http://doi.org/10.15454/1.5613785622827378E12). For multiple-label experiments, approximately 120 ng of each BAC DNA was random-priming labelled directly by incorporation of dUTP Alexa Fluor (488 or 568) or indirectly with Biotin-6-dUTP detected by immuno-FISH (Bioprime DNA labelling kit, Invitrogen, Cergy Pontoise, France). Three combinations of p or q telomeres probes of different pairs of chromosomes: (SSC2qter – SSC9qter), (SSC13qter – SSC9qter) and (SSC15qter – SSC9qter) were chosen to test their rate of association as suggested by Hi-C.
2.2.2 3D DNA Fluorescence in Situ Hybridization
3D DNA FISH experiments were conducted as described in with slight modifications. Probes were resuspended in hybridization buffer (50% formamide, 10% dextran sulphate, 2 mg/ml BSA, 2X SSC) at a final concentration of 110 ng/μL. Nuclear DNA of fibers and probes were simultaneously heat-denatured at 74°C for 7 min on the slide and then incubated overnight at 37 C in a DAKO hybridizer. Post-hybridization washes were then performed with gentle agitation, first twice in 2X SSC at 40°C for 6 min, then in 2X SSC, 50% formamide pH 7.0 at 40°C for 6 min, and finally twice for 10 min in 2X SSC, then in PBS at RT. When a biotin labelled probe was used, biotins were detected with streptavidin Alexa 568 or 488 at a final concentration of 5 μg/ml for 1 h at RT.
3D acquisitions were performed at the T.R.I. Genotoul (Toulouse Réseau Imagerie, http://trigenotoul.com/en) imaging core facility in Toulouse (France). Image stacks were collected using a Leica SP8 confocal microscope (Leica Instruments, Heidelberg, Germany) equipped with an oil immersion objective (plan achromatic 63 × N.A. = 1.4). The Z-stacks (around 80 confocal planes per capture) were acquired at 1,024 × 1,024 pixels per frame using an 8-bit pixel depth for each channel at a constant voxel size of 0.06 × 0.06 × 0.3 μm.
2.2.3 Telomere Association Analysis
Images were analyzed with specific software NEMO (), distributed under the creative commons license that can be freely downloaded from https://forge-dga.jouy.inra.fr/projects/nemo. Segmentations and 3D measurements between signals (center-to-center distance) were done as described in . Euclidean distances were computed with respect to the x, y and z resolutions. Given the resolution on the z axis, at least three pixels corresponding to 0.9 μm (0.3 × 3) were required for a high resolution between two separate signals; consequently, 1 μm was chosen as the upper cut-off for associated signals. For each combination of telomeres, nuclei were only analyzed when four signals (corresponding to the chosen telomeres of two chromosomes) were present. “Associated” signals were considered when they are separated by a distance (d) ≤ 1 μm, as done in . We determined for each combination of telomeric pairs how many nuclei were found associated among about 100 observed nuclei.
Significance of the difference in association between d90 and d110 was assessed using a χ2 test to compare generalized linear models of the binomial family with a fixed telomeric pair covariate and including, or not, the condition as a second covariate (see Supplementary Methods).
2.3 Hi-C Experiments
2.3.1 Hi-C Protocol
Hi-C experiments were performed as previously documented (), with slight modifications to adapt the Hi-C experiments and libraries to fetal muscle tissues (see Supplementary Methods).
2.3.2 Hi-C Quality Controls
After DNA digestion with HindIII, and filling-ligation of the digested ends, the HindIII target site disappears and a NheI restriction site is created instead. To check the efficiency of the Hi-C assays, PCR were performed around one HindIII restriction site with two forward primers (Fwd1: 5′ TCTGGGCAGGTCACTCATT 3′; Fwd2: 5′ TCTCGGGATGCTGAGTGTTT 3′; product size = 425 bp). A reverse primer combined with Fwd1 was used as a control (Rv1: 5′ AAACACTCAGCATCCCGAGA 3′; product size = 465 bp). In Hi-C, some religation events allow switching the sense of one DNA fragment and PCR amplification with the two forward primers. The PCR amplification products from the couple of forward primers were digested either with HindIII or NheI (product sizes = 201 + 215 bp). In control tubes (no filling of digested ends), HindIII should cleave the PCR products while NheI should not. In Hi-C tubes, NheI should cleave most of the PCR products while HindIII should cleave only a small fraction.
2.3.3 Hi-C Library Production
1.4 μg of DNA from the Hi-C experiments were fragmented with a Covaris machine. Then, 0.55 volumes of CleanPCR magnetic beads were added to the fragmented DNA to select fragments <600 bp (5 min incubation and keeping the supernatant), and 0.7 volumes of beads were added again (5 min incubation and removing supernatant) to remove fragments <200 bp. Then beads were washed with 80% ethanol and DNA was recovered with Resuspension Buffer. To purify biotinylated DNA, one volume of M-280 streptavidin magnetic Dynabeads was added and after 15 min incubation, the supernatant was removed and the beads were washed 4 times with beads wash buffer (Nextera Mate Pair Preparation Kit, Illumina) and twice with Resuspension buffer. From this point, all steps were performed while DNA remained attached to the beads. To repair DNA breaks, 60 μL of water and 40 μL of End Repair Mix 2 (TruSeqNano DNA library prep, Illumina) were added and incubated 30 min at 30°C, and then beads were washed as explained before. To allow the adapters ligation, an “A” nucleotide was added to the 3’ ends by adding 17.5 μL of water and 12.5 μL of A-Tailing Mix (TruSeqNano DNA library prep, Illumina) and incubating 30 min at 37°C and then 5 min at 70°C to inactivate the enzyme. To ligate the adapters to the DNA extremities, 2.5 μL of Resuspension Buffer, 2.5 μL of DNA Ligase Mix and 2.5 μL of DNA Adapter Index (TruSeqNano DNA library prep, Illumina) were added (10 min incubation at 30°C, then 5 μL of Stop ligation Buffer) and then beads were washed as before. DNA was amplified by 12 PCR cycles (15 s at 98°C–30 s at 60°C–30 s at 72°C) by resuspending beads in 50 μL of PCR mix (25 μL Enhanced PCR mix, 5 μL PCR primer Cocktail and 20 μL water, TruSeqNano DNA library prep, Illumina). To recover DNA from the beads, 0.6 volumes of CleanPCR magnetic beads were added and incubated 5 min, and then washed twice with 80% ethanol, resuspended in 30 μL of Resupension Buffer and after placing in a magnetic rack, supernatant containing the libraries was recovered. Libraries size was controlled with the Fragment Analyzer (FA) and quantified by qPCR. In addition, an aliquot was digested by using the NheI and HindIII enzymes to verify if selected fragments are the ones containing the filled-in biotinylated religation sites as done in . Libraries were sequenced in pool in one HiSeq3000 lane to validate their quality. For depth sequencing, the pool was paired end (PE) sequenced in 11 lines of a HiSeq3000 (reads size = 150 bases), producing from ∼ 476–685 M read pairs per library in total (see Supplementary Table 2).
2.4 Hi-C Data Analysis
2.4.1 Hi-C Reads and Interaction Matrices
The 3,447,428,742 Paired-End reads were processed using HiC-Pro v2.9.0 (
) as previously reported (
). The bioinformatics analysis includes the following steps (see
Supplementary Methodsfor more details).
• Read mapping was performed on the Sscrofa11.1 genome assembly version using Bowtie two v2.3.3.1 ().
• Interaction matrices were generated from valid pairs at various resolutions depending on the bin size. Most of the subsequent analyses were performed at the 500 Kb resolution apart from few exceptions (TAD detection for instance was performed at the 50 Kb resolution). A total of six interaction matrices were obtained per resolution (n = 3 (replicates) × 2 (groups)). Additionally, merged interaction matrices were computed by summing the interaction values of the three matrices for each group. Considering the high number of unassembled scaffolds in the pig genome Sscrofa11.1 version and given the fact that samples from both genders were collected, we focused our analysis on the 18 assembled autosomes to avoid potential effects of the sexual chromosomes on the results.
• Interaction matrices were displayed using Juicebox () and HiTC R/Bioconductor package v1.18.1 ().
• Interaction matrices were normalized per chromosome using the non-parametric iterative correction and eigenvector decomposition (ICE) method when needed ().
• Replicability between interaction matrices was assessed using the replicability index of as implemented in the R/Bioconductor package hicrep.
• Maximal resolution was computed following : a given resolution (bin size) can be claimed if, at that resolution, 80% of the bins or more contain at least 1,000 interactions. The proportion of bins with a cumulated number of valid interactions higher than 1,000 was therefore computed for different resolutions (from 100 to 5 Kb) for each individual (sample) and for the merged (group) matrix.
2.4.2 TADs Calling and Comparison
TADs were predicted per chromosome from raw interaction matrices [n = 3 (replicates) × 2 (groups) × 18 (autosomes)] at 50 Kb resolution with the Arrowhead method of the Juicer tool v1.5.3, using the -k KR parameter to ensure matrix balancing normalization. TAD finding was performed on individual matrices of each replicate separately (to assess group replicability) and on the merged matrices [n = 2 (groups) × 18 (autosomes)] to obtain a set of TADs for each group (90/110 days of gestation). To identify TADs that are consistently predicted from different replicates and group-specific TADs, we performed pairwise comparisons of TAD sets from different replicates using bedtools (v2.26.0). A mutual overlap of 90% similarity was required with the parameters -f 0.9 -r.
Insulation capacity of TAD boundaries was computed as previously described () using the local interaction score. In brief, considering all valid interactions around the same TAD boundary (i.e., both reads being not further than 500 Kb from the boundary) the interaction score corresponds to the proportion of valid interactions across the boundary. IS scores were normalized by cyclic loess () using csaw () (see Supplementary Methods for more details).
2.4.3 CTCF Prediction
The position specific frequency matrix corresponding to the CTCF-binding motif was recovered from the JASPAR Transcription Factor Binding Sites (TFBS) catalogue [http://jaspar.genereg.net, ]. CTCF genomic occurrences were predicted by running FIMO v.4.11.1 with the JASPAR CTCF frequency matrix on the Sscrofa11.1 genome. Then, the average density of CTCF predicted motifs with respect to TAD positions was obtained using bedtools v2.26.0 map and coverage functions .
2.4.4 A/B Compartments Detection
A and B compartments were obtained using the PCA approach described in , as implemented in the R/Bioconductor package HiTC (). A/B compartment identification was performed on intra-chromosome interaction matrices at 500 Kb resolution on individual interaction matrices [n = 3 (replicates) × 2 (groups) × 18 (autosomes)] and on the merged interaction matrix (n = 18 autosomes). Boundaries between A and B compartments were identified according to the sign of the first PC (eigenvector). Bins that were not assigned to any compartment due to a lack of data in some samples were not considered in subsequent integrative analyses. As an additional control, A/B compartments were also obtained by using the eigenvalue method of the Juicer tool (), which led to similar results.
The difference between the number of compartments in the two groups was assessed with a Poisson GLM: log (yijk) ∼ α cijk + βk, with yijk the number of compartments in chromosome j from sample i in group k, cijk the total number of valid interactions in chromosome j from sample i in group k, α its estimated effect on the number of compartments, and βk the estimated effect of the group on the number of compartments, which was tested for being significantly different from 0 [test with n = 2 (groups) × 3 (samples) × 18 (chromosomes) observations].
2.4.5 Detection of Differential Interactions
A differential analysis was performed to extract interactions that were significantly differentially connected between the two groups (90 and 110 days of gestation). This analysis was performed on raw count data from the 18 autosomes at the 500 Kb resolution (the differential analysis was thus performed with two groups and
n= 3 replicates in each group). A method similar to the one described in
, with some adaptations, was used to perform this task. In brief (see
Supplementary Methodsfor more details):
• Low count interactions with less 30 reads across the six samples (5 reads per sample on average) were discarded from the analysis.
• Interaction values were normalized using a non-linear normalization method based on a fast cyclic loess algorithm implemented in the R/Bioconductor package csaw ().
• Differential analysis was performed using a Generalized Linear Model (GLM) based on the Negative Binomial (NB) distribution with a group fixed effect (two-level factor: 90/110 days). The model was estimated with the implementation of the R/Bioconductor package edgeR (; ) and log ratio tests were used to assess the significativity of the group effect on each bin pair interaction. p-values were genome wide corrected using () procedure to control the False Discovery Rate (FDR).
2.4.6 Characterization of BODIs
As a single genomic bin can be involved in multiple Differential Interactions (DI) genome-wide with various logFC values, we looked for bins with a large prevalence of interactions of the same logFC sign, either mostly positive or mostly negative. A minimum ratio of 90% of DI with the same sign was required to identify “positive” or “negative” bins, possibly indicating regions that undergo a chromatin contraction or opening, respectively. Bins with a mixture of positive and negative DI were considered as undefined. Adjacent bins with the same sign (either positive, negative, or undefined) were merged into Blocks Of Differential Interactions (BODIs). This analysis was performed considering only intra-chromosomal DIs (in cis).
To assess the existence of an enrichment of large positive and negative BODIs given the relative proportions of positive and negative individual DIs, a permutation test was performed: at each permutation, logFC values were shuffled genome-wide across DIs. The same 10:1 threshold was applied to define prevalently positive and negative bins and adjacent bins of the same type were merged to identify “expected BODIs” under the null hypothesis (no specific trend of positive/negative bins to cluster consecutively). The resulting size distributions of positive, negative and undefined BODIs were compared with that of observed BODIs, and the p-value was computed, as the number of times expected BODIs were at least as frequent as the observed ones across 100 permutations for a given size and type.
The comparison of BODIs with A/B compartments was done by computing the proportion of the positive, negative and undefined BODIs that overlapped A or B compartments in terms of genomic space. The resulting block composition was therefore obtained using the bedtools coverage function on BODIs of each size and compartments of each type. As most of the compartmentalization is stable across samples, the A/B compartments obtained on the merged general matrix was used. Since A and B compartments cover roughly the same genomic space in total, no large difference should be observed between the A and B composition of positive and negative BODIs. Significance was assessed using Fisher’s exact test between the compartment type (A/B) and the BODI types (positive/negative).
2.5 Gene Expression Integrative Analysis
2.5.1 Expression Data
Expression data were obtained from a previous transcriptome study of skeletal muscle in pig during development using microarrays (). The dataset consists of 44,368 probe measurements for 17 samples (LW animals) at two different gestational stages: eight samples at 90 days and nine samples at 110 days. A precise description of the experimental design and data collection can be found in . Normalized expression data (log2 transformed) and sample information are available in NCBI (GEO accession number GSE56301). Log2 transformed expressions and log fold change (logFC) of these expression values at 90 vs 110 days (reference time point: 90 days) were used in our integrative analyses. Since the microarray was originally designed on a former version of the pig genome, probes were remapped on the Sscrofa11.1 assembly version and further filtered (see Supplementary Methods for more details).
2.5.2 Density and Expression Level of Genes in A/B Compartments
To compare the gene density in A vs B compartments, a gene density value was first computed for each compartment by dividing the number of distinct gene IDs included in the compartment (using bedtools map) by the size of the compartment. Resulting gene density distributions were then compared between A and B compartments. Normality of the gene density was tested using Shapiro-Wilk normality test and rejected for all types of compartments in both groups (p-values < 2.2e−16 overall, for n = 349 and 322 A and B compartments respectively). Wilcoxon tests were then used to assess the significance of the difference in gene density in A vs B compartments.
To compare the average gene expression in A vs B compartment, we computed for each compartment the mean expression value of its genes using bedtools map separately for the two gestational ages. Normality of the average gene expression was tested using Shapiro-Wilk normality tests and rejected for both A and B compartments (p-values = 2.58e−5 and 1.08e−3 for n = 344 and 292 A and B compartments with at least one expressed gene, respectively). Wilcoxon tests were then used to assess the significance of the difference in gene expression in A vs B compartments.
To investigate the dynamic of expression in compartment-switching regions, we considered the logFC expression values of the genes and split them into compartment-switching categories using bedtools: no switch, A to B, B to A. Normality of the logFC expression values was tested using Shapiro-Wilk normality tests for genes in all types of compartments except for compartments with no switch (n = 7, 511 genes in these compartments, above the applicability condition of the test) and rejected for both types of compartments (p-values = 1.2e−3 and 4.6e−6, for n = 60 and 174 genes in compartments switching from A to B and from B to A, respectively). Wilcoxon tests were then used to assess the significance of the difference in logFC expression values in each compartment type.
3 Results
3.1 Genome-wide Maps of Chromosomal Interactions in Fetal Porcine Muscle Tissue
We produced and sequenced Hi-C libraries from muscle samples of six pig fetuses (Supplementary Table 1): three replicates at 90 days of gestation (“d90” group) and three replicates at 110 days of gestation (“d110” group). We obtained ∼7 billion reads in total across the six samples. After trimming the sequences when needed, we could map from 63 to 73% of the pairs on the Sus scrofa v11.1 reference genome (Supplementary Table 2). These proportions are lower than usually reported with human or mouse cells (). This could be explained by several reasons, including the slightly lower quality of the porcine genomic sequence compared with the human or murine ones, and the nature of the biological material used here (frozen samples of fetal muscle). In each library, nevertheless, most of the mapped pairs showed consistent mapping configurations with respect to the genomic positions of the HindIII restriction sites (). Those were labeled as “valid interactions” (Supplementary Table 2). Overall, we obtained between 112 and 260 M valid interactions per sample from which we generated six individual interaction matrices, one per sample (Figure 1). To precisely assess the general similarity between matrices, we computed the replicability index () between all pairs of matrices from different groups (i.e., d90 vs d110) and from the same group (see Hi-C reads and interaction matrices and Supplementary Methods). By considering matrices from a previous study made on liver samples in adult pigs (), we could also compute the similarity measure between matrices from different tissues and development stages. As expected, the highest replicability index was obtained between replicates from the same group (0.92 on average, compared to 0.87 between groups and 0.67 between tissues). Adding counts from matrices of the same group generated two high-density matrices named “merged90” and “merged110” (Figure 1). More precisely, maximum matrix resolutions as defined by were 25 Kb on average per individual sample, 15 Kb for the merged110 matrix and 10 Kb for the merged90 matrix (see Methods).
FIGURE 1
3.2 A Complex Landscape of Stable and Group-specific TADs
We looked for Topologically Associating Domains (TADs) in each interaction matrix (see TADs calling and comparison) and identified 1,312 TADs per sample on average, with 84.7% of the genome being part of a TAD in at least one of the samples. Examples are displayed in Figure 2. The median TAD size of 1,200 Kb (Supplementary Table 2) was consistent with previous results in human and mouse (; Zufferey et al., 2018). In addition, computationally-predicted CTCF binding sites accumulated at TAD extremities in the expected orientation (Figure 3A, ).
FIGURE 2
FIGURE 3
The number of TADs differed between samples (from 951 to 1,585 per sample and up to 1,985 in the merged90 matrix, Supplementary Table 2). Part of this variability could be explained by the difference in the number of interactions per matrix. Indeed, computational TAD detection is known to be sensitive to variations in matrix density that can result from differences in sequencing data quantity or library complexity for instance ; Zufferey et al. (2018). Consistently, we observed a significant correlation between the number of valid interactions and the number of TADs across samples (Pearson correlation coefficient = 0.83, p-value = 9.10e−3, Figure 3B). The position of the TADs also differed between samples, although the global structure appeared stable overall (Figure 2). TAD conservation across cell types and model species has been widely reported and discussed with various degrees of circumspection (; ; ; ; ; ; ). Here, we sought to investigate TAD stability within the same tissue, by comparing samples from either the same or different groups (d90 vs d110). We considered that two TADs were identical when they overlapped with each other by 90% of their length. Pairwise comparisons of samples from the same group resulted in 1,785 identical TADs out of 2,625 on average (68.0%). As expected, this proportion was lower when comparing samples from different groups, with 1,457 identical TADs out of 2,625 on average (55.5%). Nevertheless, the observation that most of the TADs are shared within any pair of samples seemed to confirm a general stability of the TAD structure. This stability decreased drastically when requiring identity across more than two samples: for instance, only 29.0% of the TADs (2,286 out of 7,874) were identical across all the six samples. Even accounting for the presence of samples from different groups, this observed variability within the same tissue illustrates the issue of estimating TAD stability using a limited number of replicates (). The set of identical TADs in all six samples is provided in Supplementary File 1.
The difference between the proportions of identical TADs in samples from the same vs from different groups prompted us to investigate the existence of “group-specific” TADs. To find them, we considered all TADs with an identical TAD in each of the three replicates within the same group but no identical TAD in any replicate from the other group. This simple filtering process led us to a small set of 252 distinct group-specific TADs (201 for d90 and 51 for d110). It should be noted that visual inspection of the interaction matrices at the corresponding genomic positions did not show striking differences in the TAD patterns between groups (Figure 2). In order to confirm the consistency between the group-specific TADs and the raw matrix data, we computed and compared the local Interaction Score of the group-specific TAD boundaries in both groups. The Interaction Score (IS) is defined as the proportion of interactions across the midpoint of a given genomic region (see TADs calling and comparison) and can be used to assess the insulation property of TAD boundaries (). We computed the IS at each TAD boundary for each sample and computed the difference of the mean score between the d90 and the d110 groups (hereafter referred to as “ΔIS”). Negative ΔIS indicates a relative loss of interactions between 90 and 110 days. They should therefore characterize TAD boundaries that became stronger or that appeared during gestation, as one would expect for d110-specific TADs. Symmetrically, positive ΔIS indicates a gain of interactions and should therefore characterize TAD boundaries that became more permissive or disappeared during gestation. As expected, comparing the ΔIS values of the d90-and d110-specific TAD boundaries showed that the average ΔIS was positive for boundaries of d90-specific TADs but negative for boundaries of d110-specific TADs (Figure 3C). Moreover, the difference was statistically significant (p-value < 2e−7, Wilcoxon test), supporting that group-specific TADs exhibit opposite dynamics of boundary strength regardless of their number.
Considering the drastic impact TAD boundary variations can have on development (), the TAD structure differences that we observed between 90 and 110 days of gestation are likely to regulate the expression of genes involved in pig muscle maturation. Notably, we found several genes with muscle-related functions in the regions that differ between overlapping group-specific TADs, including GAP43, PECR and STIM2 for instance (; ; ). The set of group-specific TADs is provided in Supplementary File 2.
Altogether, these results showed that, while most of the TADs were preserved when comparing samples pairwise, a subset of the TADs was exclusively and consistently detected within either the d90 or the d110 group. The difference in the insulation capacity of their boundaries during time suggests that these TADs contribute to reshaping the structural organization of the pig genome during gestation.
3.3 Genome Compartments Identification Revealed a Major Chromatin Remodeling During Muscle Maturity in Pig
At a higher level of organization, we investigated the segmentation of the chromosomes into A and B epigenetic compartments using the interaction matrix of each replicate. We identified 682 compartments per replicate on average (Supplementary Table 2) with a median size between 2.6 and 3.5 Mb, in the same order of magnitude than what was reported in human or mouse cells (; ). As observed with TADs, compartment predictions were highly similar between matrices: 83.3% of the genomic regions with a prediction in each of the six samples were assigned the same compartment type in all of them consistently (either A or B six times, Figure 2 and Supplementary Figure 1), which is significantly higher than expected by chance (p-value < 1e−3, permutation test). These results illustrate the high level of reproducibility between replicates and argue for a general conservation of the higher structural organization level of the genome, as previously observed in other organisms (; ; ).
Despite this general consistency, a striking discrepancy appeared between groups. Indeed, for all replicates, d110 compartments were systematically smaller and more abundant than d90 compartments, with an increase of about 30.2% (from 593 to 772 compartments on average). A similar trend was obtained by analyzing the merged matrices (from 601 to 804 compartments for merged90 and merged110 respectively, Supplementary Table 2). This difference in the number of compartments was observed genome-wide and for both compartment types, suggesting a general fragmentation of the compartmentalization during development (Figure 4A, Supplementary Figure 1 and Supplementary Table 2). Interestingly, contrary to what was observed for TADs, no substantial correlation was detected between the total number of interactions and that of compartments (Pearson coefficient of correlation = − 0.09, p-value = 0.84, Figure 4B), ruling out variation in matrix density as a plausible explanation for this difference. These results support the idea of a major functional switch taking place in muscle cells during the maturity process, as already evidenced by expression networks () and metabolomic analyses (). Moreover, they strongly suggest that the underlying regulatory program involves epigenetic modifications through a genome-wide chromatin structure remodeling.
FIGURE 4
To investigate the potential role of such remodeling, we used gene expression data from a previous study on muscle samples at 90 and 110 days of gestation (
3.4 Comparative Analysis of Hi-C Maps Identified Significantly Different Interactions Between Gestational Stages
We then performed a comparative analysis of the Hi-C matrices to identify pairs of genomic regions with significantly different interaction values between groups of samples (see Detection of differential interactions). This analysis led to the identification of 10,183 differential interactions between pairs of 500 Kb genomic regions. While this only represents 0.11% of the 9,262,199 tested interactions, the corresponding regions involved a substantial proportion of the genomic space across all chromosomes (Figure 5A). Among the differential interactions, 8,332 (81.8%) were cis interactions, i.e., between two genomic regions from the same chromosome. This predominance is likely due to the fact that Hi-C matrices typically feature relatively low values for trans interactions, resulting in a weaker statistical power than for cis interactions.
FIGURE 5

Pairs of genomic regions with differential interactions between 90 and 110 days of gestation. Results of the comparative analysis of the Hi-C matrices at 500 Kb resolution show differential interactions along the 18 assembled autosomes (A) Differential interaction matrix. Each dot represents a pair of genomic interval with a significantly different interaction value and its associated log-fold change value (logFC, blue-white-red gradient scale). Positive values of logFC correspond to genomic regions closer at 110 days of gestation than at 90 days (red dots). Inversely, negative values indicate regions that were closer at 90 days (blue dots). Same colors are used to display cis(B) and trans(C) differential interactions as red (positive logFC) or blue (negative logFC) connections between genomic regions (outer circle). Chromosome inner color shows the genomic segmentation into A (turquoise) and B (orange) stable compartments. See also Supplementary Figures S3 and S4.
About 57% of the differential interactions showed a positive log-fold change (logFC), meaning that they contain significantly more connections at 110 days than at 90 days. These regions are therefore expected to become closer together during the 90–110 days transition. Inversely, negative logFC should characterize pairs of regions that become more distant during development. Interestingly, despite a rather balanced ratio of positive/negative logFC overall, the proportion of differential interactions with positive and negative logFC was highly heterogeneous across chromosomes (Figure 5A).
3.5 Regions Involved in Differential cis Interactions Form Homogenous Blocks of Chromatin Compaction
To further investigate the genomic distribution of significantly different interactions, we first focused on cis differential interactions and represented them along the chromosomes depending on the sign of their logFC (Figure 5B). Although each single genomic locus could potentially be involved in differential interactions of opposite logFC signs (by moving from one region to another one for instance), we noted a general predominance of one of the signs. More precisely, out of the 3,616 distinct 500 Kb regions involved in at least one differential interaction, 2,261 of them (62.5%) have either only one type (with positive or negative logFC) of interaction or at least 10 times more interactions of one type. Interestingly, regions with such a predominance of one sign tended to cluster adjacently along the genome to form homogenous blocks of either positive or negative differential interactions (Figure 5B). For instance, chromosomes 1, 13 and the q arm of chromosome 2 were largely covered by blocks of positive logFC, while blocks of negative logFC could be found in large chunks of chromosomes 3, 12 and 14 (Figure 5B). We termed these homogenous blocks BODIs, for Blocks Of Differential Interactions, and assigned to each of them its predominant logFC sign.
We first wanted to assess the significance of this observation, considering that some of the differential interactions with the same sign were expected to involve adjacent regions just by chance, necessarily forming blocks of variable sizes. To do so, we compared the size distribution of the observed BODIs with that of artificial BODIs obtained after randomly shuffling the logFC signs of the existing interactions (see Characterization of BODIs). We found a significant overrepresentation of both positive and negative BODIs of size equal or greater than 2.5 Mb up to 5 Mb (p-value < 10e−3, permutation test, Supplementary Figure 3), supporting the relevance of the observed BODIs.
Assuming that a drastic accumulation or depletion of pairwise interactions could result from variations of chromatin density, we hypothesized that positive BODIs could indicate genomic regions that undergo chromatin compaction during development. Inversely, negative BODIs would then reflect wide de-condensation events along the chromosomes. We therefore checked for consistency with the positions of A/B compartments. Interestingly, while BODIs could be found in every chromosome with a variable proportion of positive/negative BODIs, their genomic distribution in A and B compartments seemed to depend on their sign. Indeed, 58% of the genomic space in negative BODIs belonged to A compartments, while this overlap was only 30% for positive BODIs. Considering that A and B compartments occupy about the same size of the genome, this discrepancy between A and B compositions of BODIs was highly significant (p-value < 2.2e−16, Fisher’s Exact test). Consistently, a significant difference could be observed between gene expression ratios too: genes in negative BODIs had significantly higher logFC values on average than genes in positive BODIs (p-value < 2.4e−4, Wilcoxon test). These results support an epigenetic control of the chromatin compaction during late development in muscle cells.
3.6 Preferential Clustering of Telomeres at 90 Days of Gestation
We then focused on the genomic distribution of trans interactions genome-wide and observed an accumulation of differential interactions at the chromosome extremities, in particular with negative logFC (Figure 5C and Supplementary Figure 4). These interactions involved telomeric and sub-telomeric regions from both “q” and “p” arms of several chromosomes, providing additional support for a major reorganization of the chromosome conformation during gestation.
In order to validate this model, three combinations of “p” or “q” telomeric associations between different chromosomes (SSC2pter – SSC9qter, SSC13qter – SSC9qter and SSC15qter – SSC9qter) were selected based on the density of differential interactions in trans (Figures 5C and Supplementary Figure S4) and further tested by 3D DNA FISH. The number and proportion of nuclei with telomere associations were determined for each combination at 90 and 110 days. Results are presented in Figure 6 and Table 1. All three tested combinations revealed telomere clustering at both stages. Furthermore, for each combination, we obtained significantly higher proportions of association at 90 vs 110 days (p-value = 0.02, χ2 test), confirming a consistent variation of the distance between these telomeres during late gestation (Figure 6 and Table 1).
FIGURE 6

3D DNA FISH validation of preferential associations of telomeres in muscle cells. 3D images illustrating telomeric associations (SSC2pter – SSC9qter) (SSC13qter – SSC9qter) and (SSC15qter – SSC9qter) at 90 days of gestation. Maximum intensity projections of confocal image stacks are shown. SSC2p, SSC13q and SSC15q telomeres are labelled in green (Alexa 488) and SSC9qter telomere probe in red (Alexa 568). Nucleus DNA was counterstained in blue with DAPI.
TABLE 1
| Tested telomere interaction | Proportion of nuclei with interaction (total number of nuclei) | |
|---|---|---|
| 90 days of gestation | 110 days of gestation | |
| SSC9qter – SSC2pter | 24% (100) | 15% (100) |
| SSC9qter – SSC13qter | 19% (99) | 15% (100) |
| SSC9qter – SSC15qter | 28% (100) | 20% (97) |
Numbers and proportions of nuclei with an observed association between telomeres. Proportions of nuclei harboring the probed telomeric associations in muscle cells at 90 and 110 days of gestation: SSC2pter – SSC9qter, SSC13qter – SSC9qter and SSC15qter – SSC9qter. For each association, ∼100 nuclei were analyzed. A higher percentage of association is observed at 90 days of gestation for the three tested associations.
4 Discussion
4.1 First Insights in Porcine Muscle Genome Architecture During Late Gestation
To the best of our knowledge, the present study is the first 3D genome structure assessment performed on fetal muscle tissue in pig. The specific focus on the period of 14 and 4 days before birth, a critical gestation time for piglet survival at birth, makes our experimental design of high relevance for agronomic research (
As in many studies using Hi-C assays, an obvious limitation of our experimental design is the relatively low number of biological replicates, compared for instance with differential gene expression studies. Considering the ongoing cost reduction of preparing and sequencing Hi-C libraries, we expect the average number of replicates in Hi-C studies to increase in the future, as it has been the case for RNA-seq (
4.2 TAD Stability vs Variability: An Open Question
Numerous studies have led to the widespread perception that TADs are highly conserved across cell types and species (
Here, we showed that, by analyzing six samples from two different development stages of the same tissue, we could survey a wide spectrum of the topological landscape. On the one hand, pairwise comparisons between replicates of the same tissue–even from different gestational stages–resulted in a majority of identical TADs, thereby supporting the idea of a stable topological landscape. Moreover, we could identify a subset of highly stable TADs that were consistently detected in all samples regardless of the group. On the other hand, only a small proportion of the TADs (less than one third) fell into this category, meaning that the vast majority could not be found in all the samples. Also, we could identify a subset of variable TADs that were consistently group-specific, potentially enabling regulatory programs of gene expression. The presence of several genes with muscle-related function in the variable regions of these group-specific TADs supports this hypothesis, and provides interesting candidates for further functional investigations. Besides transcriptional regulation, part of this TAD variability could also be due to mechanisms like DNA replication and repair (
Overall, due to the limited relevance of any general statement on TAD variability/stability, the main challenge is probably less about estimating how variable/stable TADs are than about identifying which TADs can reliably be considered as variable/stable. In this context, ongoing efforts in data production and analysis are providing substantial help to complete and explore the known panorama of chromatin topologies, including in farm species (
4.3 Switching Compartments in Muscle Nuclei During Late Gestation
We confirmed several known features of A/B genome compartments related to gene density, expression, and general stability across replicates (
4.4 Dynamic Interacting Genomic Regions During the Maturity Process of Fetal Muscle
In this study, we could detect genome-wide dynamic changes in the chromatin structure of muscle nuclei occurring at late gestation. Specifically, we identified 10,183 differential interactions at 500 Kb resolution between the 90th and the 110th day of gestation. As noted above, considering our model of differentiated muscle fibers at two relatively close developmental stages, minor differences could have been expected. For instance, we detected much more differentially interacting regions compared with the murine myogenesis in vitro model (
The differential interactions were distributed all over the genome but not homogeneously. We observed large genomic regions of adjacent differential interactions with the same dynamic behavior when comparing the two gestational ages, sometimes along entire chromosome arms. Similar results were observed on the fly genome, where higher-order clusters corresponding to each chromosome arm were organized into active and inactive clusters (
Interestingly, we found that 58% of the genomic space in the negative BODIs was located in A compartments compared with only 30% for positive BODIs. To explain these results, we hypothesize that the structural and functional environment of A and B compartments may induce changes on the chromatin state of local regions located inside each compartment type. For instance, following our definition that negative BODIs are genomic regions that were closer (more condensed) at 90 days of gestation and that become farther apart at the end of gestation, we propose that those negative BODIs located on a decondensed/active environment (A compartment) follow a chromatin activation/de-condensation through development promoted by the genomic active environment.
4.5 Inter-chromosomal Telomeres Clustering
We found multiple dynamic associations between the telomeric regions (telomeres clustering) of several chromosomes involving either the p or the q arm. The density of interactions between telomeres decreases at 110 days of gestation. Nevertheless, 3D DNA FISH analyses do not suggest a dissociation of the clusters at the end of gestation but a higher prevalence of telomeres clustering at 90 days of gestation compared with 110 days. This indicates that telomeric regions exhibit a dynamic coordinated nuclear organization in muscle cells during late development. In fact, telomeres have been observed to display rapid movements in live human cells (
Interactions between telomeric regions have been widely reported in several species: preferential contacts between telomeres have been reported in fly embryonic nuclei, although these contacts were not associated with dynamic changes (
Interestingly, in human myoblasts, long telomeres have been observed to be involved in forming chromosome loops that can affect the higher order chromatin structure and gene expression (
In conclusion, we found major changes of the 3D genome structure during the establishment of muscle maturity at late gestation. These changes occur concomitantly with previously reported modifications of the transcriptional program, between 90 and 110 days of gestation. The topological reorganization that we observed implies structures of various scales, including individual interactions, TADs and large BODIs. The proportion of the genome that was impacted depended on the nature of the modification. Some of the changes, such as the fragmentation of the genomic A/B compartments, impacted most of the chromosomes, while others, such as the telomere clustering, involved specific regions. The amplitude of these modifications is particularly striking considering that two close fetal development stages were compared. This suggests that topological changes of the 3D genome of organized tissues could be as remarkable as changes observed during cell differentiation and cell commitment.
Statements
Data availability statement
Animals and samples metadata have been registered at BioSamples (https://www.ebi.ac.uk/biosamples) and are available using the accession ID SAMEA7390788. Experimental protocols have been deposited at the FAANG DCC: https://data.faang.org/api/fire_api/samples/INRAE_SOP_pig_muscle_tissue_sampling_20200812.pdf (sampling) and https://data.faang.org/api/fire_api/assays/INRAE_SOP_Hi-C_pig_muscle_tissue_20200812.pdf (Hi-C libraries). Raw sequencing data from Hi-C experiments have been uploaded to the ENA under the accession PRJEB40576 (ERP124229) and to the FAANG DCC along with the associated metadata, as per the FAANG consortium agreement (https://www.faang.org/data-share-principle). Hi-C matrices, TADs and A/B compartments generated in this study have been made available at: https://doi.org/10.15454/DOMEHB and http://www.fragencode.org/pig3Dgenome.html.
Ethics statement
The animal study was reviewed and approved by the ethical committee (#84) in animal experimentation of the French Ministry of National Education, Higher Education, and Scientific Research (authorization #02015021016014354). The experiment authorization number for the experimental farm GenESI (Genetics, testing and innovative systems experimental unit) is A17661. The procedures performed in this study and the treatment of animals complied with European Union legislation (Directive 2010/63/EU) and French legislation in the Midi-Pyrénées Region of France (Decree 2001–464).
Author contributions
MY-B, MM-M, YL-M, and SF conceived the study. MM-M and YL-M designed and supervised the experiments. MM-M performed the Hi-C experiments, YL-M and SC performed the 3D DNA FISH experiments. SF and NV designed and supervised the data analysis. SF, NV, MZ, MM-M, and DR performed the bioinformatics and statistical analyses. MM-M, SF, NV, MY-B, and YL-M wrote the manuscript with the input from all authors.
Funding
This work was supported by a grant to MM-M from the French “Ministère de l’Education Nationale, de la Recherche et de la Technologie”. Tissue sampling and Hi-C library sequencing were funded by the 2014 Animal Genetics Division AAP grant (INRA, France), and by the “Osez l’interdisciplinarité !” grant of P. Neuvial's SCALES project (CNRS, France). Animals were raised in the experimental farm GenESI (INRAE, France). Further support was provided to MM-M with a training program founded by the “Soutien à la mobilité internationale” scholarship from the INP (Toulouse, France) and hosted by M. Groenen’s group (WUR, Netherlands).
Acknowledgments
We would like to thank the GeT-PlaGe facility (INRAE, France), the T.R.I. Genotoul imaging core facility (Toulouse, France), and the Genotoul bioinformatics platform (Toulouse, France) for providing assistance and technical support. We also thank L. Liaubet’s team (INRAE, France) for providing samples and assisting with the experiments, H. Acloque (INRAE, France) for pioneering the Hi-C in our institute, and both the Chrocogen working group and the FAANG community for their support.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.748239/full#supplementary-material
References
1
BallmanK. V.GrillD. E.ObergA. L.TherneauT. M. (2004). Faster Cyclic Loess: Normalizing RNA Arrays via Linear Models. Bioinformatics20, 2778–2786. 10.1093/bioinformatics/bth327
2
BarutcuA. R.LajoieB. R.FritzA. J.McCordR. P.NickersonJ. A.Van WijnenA. J.et al (2016). SMARCA4 Regulates Gene Expression and Higher-Order Chromatin Structure in Proliferating Mammary Epithelial Cells. Genome Res.26, 1188–1201. 10.1101/gr.201624.115
3
BarutcuA. R.LajoieB. R.McCordR. P.TyeC. E.HongD.MessierT. L.et al (2015). Chromatin Interaction Analysis Reveals Changes in Small Chromosome and Telomere Clustering between Epithelial and Breast Cancer Cells. Genome Biol.16, 214. 10.1186/s13059-015-0768-0
4
BeltonJ.-M.McCordR. P.GibcusJ. H.NaumovaN.ZhanY.DekkerJ. (2012). Hi-C: a Comprehensive Technique to Capture the Conformation of Genomes. Methods58, 268–276. 10.1016/j.ymeth.2012.05.001
5
BenjaminiY.HochbergY. (1995). Controlling the False Discovery Rate: a Practical and Powerful Approach to Multiple Testing. J. R. Stat. Soc. Ser. B (Methodological)57, 289–300. 10.1111/j.2517-6161.1995.tb02031.x
6
BolzerA.KrethG.SoloveiI.KoehlerD.SaracogluK.FauthC.et al (2005). Three-dimensional Maps of All Chromosomes in Human Male Fibroblast Nuclei and Prometaphase Rosettes. Plos Biol.3, e157. 10.1371/journal.pbio.0030157
7
BonevB.CavalliG. (2016). Organization and Function of the 3D Genome. Nat. Rev. Genet.17, 661–678. 10.1038/nrg.2016.112
8
BonevB.Mendelson CohenN.SzaboQ.FritschL.PapadopoulosG. L.LublingY.et al (2017). Multiscale 3D Genome Rewiring during Mouse Neural Development. Cell171, 557–572. e24. 10.1016/j.cell.2017.09.043
9
BonoraG.PlathK.DenholtzM. (2014). A Mechanistic Link between Gene Regulation and Genome Architecture in Mammalian Development. Curr. Opin. Genet. Development27, 92–101. 10.1016/j.gde.2014.05.002
10
BoyaR.YadavalliA. D.NikhatS.KurukutiS.PalakodetiD.PongubalaJ. M. R. (2017). Developmentally Regulated Higher-Order Chromatin Interactions Orchestrate B Cell Fate Commitment. Nucleic Acids Res.45, 11070–11087. 10.1093/nar/gkx722
11
BrancoM. R.PomboA. (2006). Intermingling of Chromosome Territories in Interphase Suggests Role in Translocations and Transcription-dependent Associations. Plos Biol.4, e138. 10.1371/journal.pbio.0040138
12
CremerT.CremerC. (2001). Chromosome Territories, Nuclear Architecture and Gene Regulation in Mammalian Cells. Nat. Rev. Genet.2, 292–301. 10.1038/35066075
13
CremerT.CremerM.DietzelS.MüllerS.SoloveiI.FakanS. (2006). Chromosome Territories - a Functional Nuclear Landscape. Curr. Opin. Cell Biol.18, 307–316. 10.1016/j.ceb.2006.04.007
14
Cubeñas-PottsC.CorcesV. G. (2015). Topologically Associating Domains: an Invariant Framework or a Dynamic Scaffold. Nucleus6, 430–434. 10.1080/19491034.2015.1096467
15
DaliR.BlanchetteM. (2017). A Critical Assessment of Topologically Associating Domain Prediction Tools. Nucleic Acid Res.45, 2994–3005. 10.1093/nar/gkx145
16
DarbellayB.ArnaudeauS.CeroniD.BaderC. R.KonigS.BernheimL. (2010). Human Muscle Economy Myoblast Differentiation and Excitation-Contraction Coupling Use the Same Molecular Partners, Stim1 and Stim2. J. Biol. Chem.285, 22437–22447. 10.1074/jbc.M110.118984
17
DaviesJ. O. J.OudelaarA. M.HiggsD. R.HughesJ. R. (2017). How Best to Identify Chromosomal Interactions: a Comparison of Approaches. Nat. Methods14, 125–134. 10.1038/nmeth.4146
18
DekkerJ.RippeK.DekkerM.KlecknerN. (2002). Capturing Chromosome Conformation. Science295, 1306–1311. 10.1126/science.1067799
19
DixonJ. R.JungI.SelvarajS.ShenY.Antosiewicz-BourgetJ. E.LeeA. Y.et al (2015). Chromatin Architecture Reorganization during Stem Cell Differentiation. Nature518, 331–336. 10.1038/nature14222
20
DixonJ. R.SelvarajS.YueF.KimA.LiY.ShenY.et al (2012). Topological Domains in Mammalian Genomes Identified by Analysis of Chromatin Interactions. Nature485, 376–380. 10.1038/nature11082
21
DoynovaM. D.MarkworthJ. F.Cameron-SmithD.VickersM. H.O’SullivanJ. M. (2017). Linkages between Changes in the 3D Organization of the Genome and Transcription during Myotube Differentiation In Vitro. Skeletal Muscle7, 1–14. 10.1186/s13395-017-0122-1
22
DurandN. C.ShamimM. S.MacholI.RaoS. S. P.HuntleyM. H.LanderE. S.et al (2016). Juicer Provides a One-Click System for Analyzing Loop-Resolution Hi-C Experiments. Cell Syst.3, 95–98. 10.1016/j.cels.2016.07.002
23
EresI. E.GiladY. (2021). A TAD Skeptic: Is 3D Genome Topology Conserved. Trends Genet.37, 216–223. 10.1016/j.tig.2020.10.009
24
FeodorovaY.FalkM.MirnyL. A.SoloveiI. (2020). Viewing Nuclear Architecture through the Eyes of Nocturnal Mammals. Trends Cell Biol.30, 276–289. 10.1016/j.tcb.2019.12.008
25
FinnE. H.PegoraroG.BrandãoH. B.ValtonA.-L.OomenM. E.DekkerJ.et al (2019). Extensive Heterogeneity and Intrinsic Variation in Spatial Genome Organization. Cell176, 1502–1515. 10.1016/j.cell.2019.01.020
26
FoissacS.DjebaliS.MunyardK.VialaneixN.RauA.MuretK.et al (2019). Multi-species Annotation of Transcriptome and Chromatin Structure in Domesticated Animals. BMC Biol.17, 108. 10.1186/s12915-019-0726-5
27
FoxcroftG. R.DixonW. T.NovakS.PutmanC. T.TownS. C.VinskyM. D. A. (2006). The Biological Basis for Prenatal Programming of Postnatal Performance in Pigs1,2. J. Anim. Sci.84, E105–E112. 10.2527/2006.8413_supplE105x
28
FraserJ.FerraiC.ChiarielloA. M.SchuelerM.RitoT.LaudannoG.et al (2015). Hierarchical Folding and Reorganization of Chromosomes Are Linked to Transcriptional Changes in Cellular Differentiation. Mol. Syst. Biol.11, 852. 10.15252/msb.20156492
29
GibcusJ. H.DekkerJ. (2013). The Hierarchy of the 3D Genome. Mol. Cell49, 773–782. 10.1016/j.molcel.2013.02.011
30
GiuffraE.TuggleC. K.ConsortiumF. A. A. N. G. (2019). Functional Annotation of Animal Genomes (FAANG): Current Achievements and Roadmap. Annu. Rev. Anim. Biosci.7, 65–88. 10.1146/annurev-animal-020518-114913
31
GrantC. E.BaileyT. L.NobleW. S. (2011). FIMO: Scanning for Occurrences of a Given Motif. Bioinformatics27, 1017–1018. 10.1093/bioinformatics/btr064
32
GuarnieriS.MorabitoC.PaoliniC.BoncompagniS.PillaR.Fanò-IllicG.et al (2013). Growth Associated Protein 43 Is Expressed in Skeletal Muscle Fibers and Is Localized in Proximity of Mitochondria and Calcium Release Units. PLoS One8, e53267. 10.1371/journal.pone.0053267
33
GuidiM.RuaultM.MarboutyM.LoïodiceI.CournacA.BillaudeauC.et al (2015). Spatial Reorganization of Telomeres in Long-Lived Quiescent Cells. Genome Biol.16, 206. 10.1186/s13059-015-0766-2
34
HaradaA.MallappaC.OkadaS.ButlerJ. T.BakerS. P.LawrenceJ. B.et al (2015). Spatial Re-organization of Myogenic Regulatory Sequences Temporally Controls Gene Expression. Nucleic Acids Res.43, 2008–2021. 10.1093/nar/gkv046
35
HarewoodL.KishoreK.EldridgeM. D.WingettS.PearsonD.SchoenfelderS.et al (2017). Hi-C as a Tool for Precise Detection and Characterisation of Chromosomal Rearrangements and Copy Number Variation in Human Tumours. Genome Biol.18, 125. 10.1186/s13059-017-1253-8
36
HeM.LiY.TangQ.LiD.JinL.TianS.et al (2018). Genome-wide Chromatin Structure Changes during Adipogenesis and Myogenesis. Int. J. Biol. Sci.14, 1571–1585. 10.7150/ijbs.25328
37
HozéN.RuaultM.AmorusoC.TaddeiA.HolcmanD. (2013). Spatial Telomere Organization and Clustering in Yeast Saccharomyces cerevisiae Nucleus Is Generated by a Random Dynamics of Aggregation-Dissociation. MBoC24, 1791–1800. 10.1091/mbc.E13-01-0031
38
IannuccelliE.MompartF.GellinJ.Lahbib-MansaisY.YerleM.BoudierT. (2010). NEMO: a Tool for Analyzing Gene and Chromosome Territory Distributions from 3D-FISH Experiments. Bioinformatics26, 696–697. 10.1093/bioinformatics/btq013
39
ImakaevM.FudenbergG.McCordR. P.NaumovaN.GoloborodkoA.LajoieB. R.et al (2012). Iterative Correction of Hi-C Data Reveals Hallmarks of Chromosome Organization. Nat. Methods9, 999–1003. 10.1038/nmeth.2148
40
Lahbib-MansaisY.BarascH.Marti-MarimonM.MompartF.IannuccelliE.RobelinD.et al (2016). Expressed Alleles of Imprinted IGF2, DLK1 and MEG3 Colocalize in 3D-Preserved Nuclei of Porcine Fetal Cells. BMC Cell Biol.17, 1–15. 10.1186/s12860-016-0113-9
41
LangmeadB.SalzbergS. L. (2012). Fast Gapped-Read Alignment with Bowtie 2. Nat. Methods9, 357–359. 10.1038/nmeth.1923
42
Lazar-StefanitaL.ScolariV. F.MercyG.MullerH.GuérinT. M.ThierryA.et al (2017). Cohesins and Condensins Orchestrate the 4D Dynamics of Yeast Chromosomes during the Cell Cycle. EMBO J.36, 2684–2697. 10.15252/embj.201797342
43
LefortG.ServienR.QuesnelH.BillonY.CanarioL.IannuccelliN.et al (2020). The Maturity in Fetal Pigs Using a Multi-Fluid Metabolomic Approach. Sci. Rep.10, 19912. 10.1038/s41598-020-76709-8
44
Lieberman-AidenE.Van BerkumN. L.WilliamsL.ImakaevM.RagoczyT.TellingA.et al (2009). Comprehensive Mapping of Long-Range Interactions Reveals Folding Principles of the Human Genome. Science326, 289–293. 10.1126/science.1181369
45
LinY. C.BennerC.ManssonR.HeinzS.MiyazakiK.MiyazakiM.et al (2012). Global Changes in the Nuclear Positioning of Genes and Intra- and Interdomain Genomic Interactions that Orchestrate B Cell Fate. Nat. Immunol.13, 1196–1204. 10.1038/ni.2432
46
LiuY.ZhouJ.WhiteK. P. (2014). RNA-seq Differential Expression Studies: More Sequence or More Replication. Bioinformatics30, 301–304. 10.1093/bioinformatics/btt688
47
LongH. K.PrescottS. L.WysockaJ. (2016). Ever-changing Landscapes: Transcriptional Enhancers in Development and Evolution. Cell167, 1170–1187. 10.1016/j.cell.2016.09.018
48
LunA. T. L.SmythG. K. (2016). Csaw: a Bioconductor Package for Differential Binding Analysis of ChIP-Seq Data Using Sliding Windows. Nucleic Acids Res.44, e45. 10.1093/nar/gkv1191
49
LunA. T. L.SmythG. K. (2015). diffHic: a Bioconductor Package to Detect Differential Genomic Interactions in Hi-C Data. BMC Bioinformatics16, 258. 10.1186/s12859-015-0683-0
50
LunneyJ. K. (2007). Advances in Swine Biomedical Model Genomics. Int. J. Biol. Sci.3, 179–184. 10.7150/ijbs.3.179
51
LupiáñezD. G.KraftK.HeinrichV.KrawitzP.BrancatiF.KlopockiE.et al (2015). Disruptions of Topological Chromatin Domains Cause Pathogenic Rewiring of Gene-Enhancer Interactions. Cell161, 1012–1025. 10.1016/j.cell.2015.04.004
52
Marti-MarimonM.VialaneixN.VoilletV.Yerle-BouissouM.Lahbib-MansaisY.LiaubetL. (2018). A New Approach of Gene Co-expression Network Inference Reveals Significant Biological Processes Involved in Porcine Muscle Development in Late Gestation. Sci. Rep.8, 10150. 10.1038/s41598-018-28173-8
53
MathelierA.FornesO.ArenillasD. J.ChenC.-y.DenayG.LeeJ.et al (2016). JASPAR 2016: a Major Expansion and Update of the Open-Access Database of Transcription Factor Binding Profiles. Nucleic Acids Res.44, D110–D115. 10.1093/nar/gkv1176
54
McCarthyD. J.ChenY.SmythG. K. (2012). Differential Expression Analysis of Multifactor RNA-Seq Experiments with Respect to Biological Variation. Nucleic Acids Res.40, 4288–4297. 10.1093/nar/gks042
55
McCordR. P.KaplanN.GiorgettiL. (2020). Chromosome Conformation Capture and beyond: toward an Integrative View of Chromosome Structure and Function. Mol. Cell77, 688–708. 10.1016/j.molcel.2019.12.021
56
MeurensF.SummerfieldA.NauwynckH.SaifL.GerdtsV. (2012). The Pig: a Model for Human Infectious Diseases. Trends Microbiol.20, 50–57. 10.1016/j.tim.2011.11.002
57
MolenaarC.WiesmeijerK.VerwoerdN. P.KhazenS.EilsR.TankeH. J.et al (2003). Visualizing Telomere Dynamics in Living Mammalian Cells Using PNA Probes. EMBO J.22, 6631–6641. 10.1093/emboj/cdg633
58
MompartF.RobelinD.DelcrosC.Yerle-BouissouM. (2013). 3d Organization of Telomeres in Porcine Neutrophils and Analysis of LPS-Activation Effect. BMC Cell Biol.14, 30. 10.1186/1471-2121-14-30
59
NaganoT.LublingY.StevensT. J.SchoenfelderS.YaffeE.DeanW.et al (2013). Single-cell Hi-C Reveals Cell-To-Cell Variability in Chromosome Structure. Nature502, 59–64. 10.1038/nature12593
60
NageleR. G.VelascoA. Q.AndersonW. J.McMahonD. J.ThomsonZ.FazekasJ.et al (2001). Telomere Associations in Interphase Nuclei: Possible Role in Maintenance of Interphase Chromosome Topology. J. Cell Sci.114, 377–388. 10.1242/jcs.114.2.377
61
NoraE. P.LajoieB. R.SchulzE. G.GiorgettiL.OkamotoI.ServantN.et al (2012). Spatial Partitioning of the Regulatory Landscape of the X-Inactivation centre. Nature485, 381–385. 10.1038/nature11049
62
PiórkowskaK.TyraM.Ropka-MolikK.PodbielskaA. (2017). Evolution of Peroxisomal Trans-2-enoyl-CoA Reductase ( PECR ) as Candidate Gene for Meat Quality. Livestock Sci.201, 85–91. 10.1016/j.livsci.2017.05.004
63
QuinlanA. R. (2014). BEDTools: the Swiss-army Tool for Genome Feature Analysis. Curr. Protoc. Bioinformatics47, 11–34. 1–34. 10.1002/0471250953.bi1112s47
64
RamaniV.DengX.QiuR.GundersonK. L.SteemersF. J.DistecheC. M.et al (2017). Massively Multiplex Single-Cell Hi-C. Nat. Methods14, 263–266. 10.1038/nmeth.4155
65
RaoS. S. P.HuntleyM. H.DurandN. C.StamenovaE. K.BochkovI. D.RobinsonJ. T.et al (2014). A 3D Map of the Human Genome at Kilobase Resolution Reveals Principles of Chromatin Looping. Cell159, 1665–1680. 10.1016/j.cell.2014.11.021
66
RapaportF.KhaninR.LiangY.PirunM.KrekA.ZumboP.et al (2013). Comprehensive Evaluation of Differential Gene Expression Analysis Methods for RNA-Seq Data. Genome Biol.14, R95–R13. 10.1186/gb-2013-14-9-r95
67
RehfeldtC.FiedlerI.DietlG.EnderK. (2000). Myogenesis and Postnatal Skeletal Muscle Cell Growth as Influenced by Selection. Livestock Prod. Sci.66, 177–188. 10.1016/S0301-6226(00)00225-6
68
RehfeldtC.KuhnG. (2006). Consequences of Birth Weight for Postnatal Growth Performance and Carcass Quality in Pigs as Related to Myogenesis1. J. Anim. Sci.84, E113–E123. 10.2527/2006.8413_supplE113x
69
RobinJ. D.LudlowA. T.BattenK.MagdinierF.StadlerG.WagnerK. R.et al (2014). Telomere Position Effect: Regulation of Gene Expression with Progressive Telomere Shortening over Long Distances. Genes Dev.28, 2464–2476. 10.1101/gad.251041.114
70
RobinsonM. D.McCarthyD. J.SmythG. K. (2010). edgeR: a Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics26, 139–140. 10.1093/bioinformatics/btp616
71
SauerwaldN.SinghalA.KingsfordC. (2020). Analysis of the Structural Variability of Topologically Associated Domains as Revealed by Hi-C. NAR Genomics and Bioinformatics2, lqz008. 10.1093/nargab/lqz008
72
SchmittA. D.HuM.JungI.XuZ.QiuY.TanC. L.et al (2016). A Compendium of Chromatin Contact Maps Reveals Spatially Active Regions in the Human Genome. Cell Rep.17, 2042–2059. 10.1016/j.celrep.2016.10.061
73
ServantN.LajoieB. R.NoraE. P.GiorgettiL.ChenC.-J.HeardE.et al (2012). HiTC: Exploration of High-Throughput ′C′ Experiments. Bioinformatics28, 2843–2844. 10.1093/bioinformatics/bts521
74
ServantN.VaroquauxN.LajoieB. R.ViaraE.ChenC.-J.VertJ.-P.et al (2015). HiC-Pro: an Optimized and Flexible Pipeline for Hi-C Data Processing. Genome Biol.16, 259. 10.1186/s13059-015-0831-x
75
SextonT.YaffeE.KenigsbergE.BantigniesF.LeblancB.HoichmanM.et al (2012). Three-dimensional Folding and Functional Organization Principles of the Drosophila Genome. Cell148, 458–472. 10.1016/j.cell.2012.01.010
76
Solov’evaL.SvetlovaM.BodinskiD.ZalenskyA. O. (2004). Nature of Telomere Dimers and Chromosome Looping in Human Spermatozoa. Chromosome Res.12, 817–823. 10.1007/s10577-005-5513-1
77
StevensT. J.LandoD.BasuS.AtkinsonL. P.CaoY.LeeS. F.et al (2017). 3D Structures of Individual Mammalian Genomes Studied by Single-Cell Hi-C. Nature544, 59–64. 10.1038/nature21429
78
TanL.XingD.ChangC.-H.LiH.XieX. S. (2018). Three-dimensional Genome Structures of Single Diploid Human Cells. Science361, 924–928. 10.1126/science.aat5641
79
Vietri RudanM.BarringtonC.HendersonS.ErnstC.OdomD. T.TanayA.et al (2015). Comparative Hi-C Reveals that CTCF Underlies Evolution of Chromosomal Domain Architecture. Cell Rep10, 1297–1309. 10.1016/j.celrep.2015.02.004
80
Vietri RudanM.BarringtonC.HendersonS.ErnstC.OdomD. T.TanayA.et al (2015). Comparative Hi-C Reveals that CTCF Underlies Evolution of Chromosomal Domain Architecture. Cell Rep.10, 1297–1309. 10.1016/j.celrep.2015.02.004
81
VoilletV.San CristobalM.PèreM.-C.BillonY.CanarioL.LiaubetL.et al (2018). Integrated Analysis of Proteomic and Transcriptomic Data Highlights Late Fetal Muscle Maturation Process. Mol. Cell Proteomics17, 672–693. 10.1074/mcp.M116.066357
82
VoilletV.SanCristobalM.LippiY.MartinP. G.IannuccelliN.LascorC.et al (2014). Muscle Transcriptomic Investigation of Late Fetal Development Identifies Candidate Genes for Piglet Maturity. BMC Genomics15, 797. 10.1186/1471-2164-15-797
83
WangX.KamZ.CarltonP. M.XuL.SedatJ. W.BlackburnE. H. (2008). Rapid Telomere Motions in Live Human Cells Analyzed by Highly Time-Resolved Microscopy. Epigenetics & Chromatin1, 4. 10.1186/1756-8935-1-4
84
WonH.de La Torre-UbietaL.SteinJ. L.ParikshakN. N.HuangJ.OplandC. K.et al (2016). Chromosome Conformation Elucidates Regulatory Relationships in Developing Human Brain. Nature538, 523–527. 10.1038/nature19847
85
YaffeE.TanayA. (2011). Probabilistic Modeling of Hi-C Contact Maps Eliminates Systematic Biases to Characterize Global Chromosomal Architecture. Nat. Genet.43, 1059–1065. 10.1038/ng.947
86
YamamotoA. (2014). Gathering up Meiotic Telomeres: a Novel Function of the Microtubule-Organizing center. Cell. Mol. Life Sci.71, 2119–2134. 10.1007/s00018-013-1548-1
87
YangT.ZhangF.YardımcıG. G.SongF.HardisonR. C.NobleW. S.et al (2017). HiCRep: Assessing the Reproducibility of Hi-C Data Using a Stratum-Adjusted Correlation Coefficient. Genome Res.27, 1939–1949. 10.1101/gr.220640.117
88
Yerle-BouissouM.MompartF.IannuccelliE.RobelinD.JauneauA.Lahbib-MansaisY.et al (2009). Nuclear Architecture of Resting and LPS-Stimulated Porcine Neutrophils by 3D FISH. Chromosome Res.17, 847–862. 10.1007/s10577-009-9074-6
89
ZhangN.Mendieta-EstebanJ.MagliA.LiljaK. C.PerlingeiroR. C. R.Marti-RenomM. A.et al (2020). Muscle Progenitor Specification and Myogenic Differentiation Are Associated with Changes in Chromatin Topology. Nat. Commun.11, 6222. 10.1038/s41467-020-19999-w
90
ZhengH.XieW. (2019). The Role of 3d Genome Organization in Development and Cell Differentiation. Nat. Rev. Mol. Cell Biol.20, 535–550. 10.1038/s41580-019-0132-4
91
ZuffereyM.TavernariD.OricchioE.CirielloG. (2018). Comparison of Computational Methods for the Identification of Topologically Associating Domains. Genome Biol.19, 217. 10.1186/s13059-018-1596-9
Summary
Keywords
Hi-C, chromosome conformation, chromatin structure, telomeres, pig, development, fetal muscle, 3D genome architecture
Citation
Marti-Marimon M, Vialaneix N, Lahbib-Mansais Y, Zytnicki M, Camut S, Robelin D, Yerle-Bouissou M and Foissac S (2021) Major Reorganization of Chromosome Conformation During Muscle Development in Pig. Front. Genet. 12:748239. doi: 10.3389/fgene.2021.748239
Received
27 July 2021
Accepted
14 September 2021
Published
05 October 2021
Volume
12 - 2021
Edited by
Francisco Peñagaricano, University of Wisconsin-Madison, United States
Reviewed by
Cong-jun Li, Agricultural Research Service (USDA), United States
Tosso Leeb, University of Bern, Switzerland
Wai Yee Low, University of Adelaide, Australia
Updates

Check for updates
Copyright
© 2021 Marti-Marimon, Vialaneix, Lahbib-Mansais, Zytnicki, Camut, Robelin, Yerle-Bouissou and Foissac.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Sylvain Foissac, sylvain.foissac@inrae.fr
† Present address: CNAG-CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), Baldiri i Reixac 4, 08,028 Barcelona, Spain
This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.