Dynamic transcriptome and chromatin architecture in granulosa cells during chicken folliculogenesis

Folliculogenesis is a complex biological process involving a central oocyte and its surrounding somatic cells. Three-dimensional chromatin architecture is an important transcription regulator; however, little is known about its dynamics and role in transcriptional regulation of granulosa cells during chicken folliculogenesis. We investigate the transcriptomic dynamics of chicken granulosa cells over ten follicular stages and assess the chromatin architecture dynamics and how it influences gene expression in granulosa cells at three key stages: the prehierarchical small white follicles, the first largest preovulatory follicles, and the postovulatory follicles. Our results demonstrate the consistency between the global reprogramming of chromatin architecture and the transcriptomic divergence during folliculogenesis, providing ample evidence for compartmentalization rearrangement, variable organization of topologically associating domains, and rewiring of the long-range interaction between promoter and enhancers. These results provide key insights into avian reproductive biology and provide a foundational dataset for the future in-depth functional characterization of granulosa cells.

compartment, TAD and PEI) respectively contribute to changes in gene expression. 2. RNA-seq readout corresponds to mRNA levels. Both transcriptional activities and RNA degradation affect mRNA levels. The authors should discuss the relationships between changes in chromatin compartmentalisation and alterations in mRNA levels. 3. The authors should add a two-dimensional scatter plot to visualise the global relationships among RNA-seq samples 4. The authors have used maSigPro to cluster all genes to globally characterise transcriptomic changes and then used GO analysis to summarise the functions of genes assigned to each cluster. A well-known limitation of GO analysis is its over-generalisation and vagueness. Such an example is, for instance, the enrichment in the "growth and developmental processes" GO term. The authors should include a heatmap of enriched/representative genes that define particular terms between stages to replace fig 1d. Another option would be to use bubble plots to include both GO term enrichment and statistical significance values in the exact representation. 5. Given the study's descriptive nature and lack of resolution in chromatin data, the authors should define specific active putative cis-regulatory elements by adding chromatin accessibility (ATACseq) characterisation at selected time points (same as Hi-C analysis). Integration of these data would greatly help understand the intricacies of the gene regulatory programme underlying chicken ovarian follicle development. 6. It is surprising that using pairwise differential expression analysis, the authors found "no difference in gene expression between contiguous stages of development" (line 139-140), especially between SWF/LWF, LYF/F5 and F1/POF. Increasing the sequencing depth would likely for RNA-seq libraries as well allow identification of the differentially expressed genes between these stages (see comment above) 7. As admitted by the authors, "chromatin compartmentalisation contributes to relatively subtle changes in gene expression and does not play a deterministic role" (line 314-317). Hence, by definition, the significance of this work in terms of understanding the determinants of gene expression during follicle development is limited. One way to possibly improve the study's depth, bringing new insights, would be by extending the transcriptomic analysis using single-cell RNAseq on selected developmental stages -the study would then provide an atlas of both transcriptome and genome-wide chromatin interactions during chicken ovarian follicle development.

Detailed responses to reviewers
All comments provided by reviewers are given in gray italics, and our responses are in black. All revisions in the manuscript are marked in red.

Reviewer #1:
In the manuscript "Dynamic chromatin architecture in ovarian follicle development in chickens" authors provided a detailed analysis of gene We took steps to control the purity of granulosa cells (GCs). First, the follicular GCs were collected according to the methods described by Gilbert. et al. (Gilbert et al., 1977). We have added the below descriptions to the Methods section of the main text (Main text page 39, lines 650-661). The follicles were carefully excised from the ovaries of birds under general anesthesia or euthanized with sodium pentobarbitone. The stalk of the excised follicle was held with a forceps so that the clear, avascular stigma was on top and a cut was made with a scalpel approximately along the line of the stigma quickly.
Immediately after it was cut, the follicle was inverted over phosphate buffer solution (PBS) and the follicles were shaken to remove the yolk. The follicles were gently shaken until a transparent film appeared in the PBS, which is the granular cell.
Furthermore, to demonstrate the purity of the GC samples, we constructed single-cell libraries for GCs at the SWF, F1, and POF stages. After quality control and filtering, a total of 21,393 high-quality cells (6,596, 5,996, and 8,801 cells for SWF, F1, and POF stages, respectively) were available for cell-type characterization. By integrating scRNA-seq data of developing chicken hearts with our data, we classified the cells into 13 clusters using the unbiased clustering and uniform manifold approximation and projection (UMAP) methods, generated hierarchical clustering using the 50 most variably expressed gene means per cluster, and finally distinguished two major cell groups: the GC group (three clusters) and other cells (ten clusters). To confirm the identity of these GCs, we colored the single cells according to the expression levels of five canonical markers of GCs (CYP11A1, CHST8, FSHR, TSPAN6, and DSP) and five representative genes (NOV, RLN3, EDN2, FGL2, and RGS16) specifically high expressed in GCs. We found that almost all (>99%) sequenced cells at the three stages were annotated as GCs. These results suggested that the isolated GCs had little contamination and high purity.

Comment 1-2:
"Focusing on the TAD boundaries across the development of follicles, the intensified spatial segregation at the F1 stage indicated that more selfinteractions occurred within TADs when the follicle reached maturity (Fig. 4df)." -the portion of intra-TAD interactions may reflect noise level and/or P(s) scaling, which is slightly different between samples. For example, it is expectable that higher noise level will lead to higher portion of intra-TAD interactions. This possibility should be excluded before drawing biological conclusions from very slight differences shown in Fig. 4 d and e. For example, note that in Fig. 4,E at the point -300 kb, the difference in insulation between POF and other stages is higher than in the TAD boundary region.

Response 1-2:
Thank you very much for these helpful comments. To effectively reduce the noise compared to intra-TAD interactions, we normalized TAD lengths by calculating TAD interaction enrichment and TAD intactness (Figure 4 e,f) in our revised manuscript. After normalization, the spatial segregation at the F1 stage is significantly higher than at the SWF stage (P < 0.001) and the POF stage (P < 2.2×10 -16 ). Furthermore, the entropy status (VNE) of a cis-contact matrix was explored and confirms this result (Figure 2a As for the previous Figure 4e version, it seems normal that the difference in insulation between the POF stage and other stages is higher than in the TAD boundary region at the point -300 kb. When identifying the TAD boundary according to the insulation score (IS) algorithm, we normalized the IS of a bin relative to all bins across that chromosome using the following formula: log2(IS/mean_IS(all_bins)). Valleys/minima along the normalized IS vector represent loci of reduced Hi-C interactions across the bins and are considered TAD boundaries. As such, in our study, the median length of the TAD boundary is only 20 kb (at 20 kb resolution). The difference at the point ±300 kb away from the center of the TAD boundary does not reflect the main characteristics of TADs, nor does it conflict with our conclusion about intra-TAD interactions.

Comment 1-3：
Promoter-enhancer interactions should be distinguished from CTCF-based loops formed via loop extrusion mechanism. It is known that CTCF sites are preferentially located near promoters and/or enhancers. Consistent with this, it seems that most of interactions depicted in Fig. 6, G and called by PSYCHIC connect TAD boundaries. These CTCF-mediated interactions might have regulatory effect; however, acute degradation of CTCF protein (Rao et al., 2017) showed that there are "bona fide" PEIs, which are mediated by mechanisms other than CTCF blocked extrusion. Authors should carefully distinguish these two types of interactions through the whole paragraph "Global remodeling of promoter-enhancer interactions in transcription control during follicle development".

Response 1-3:
Thank you for this suggestion. We tried to construct CTCF ChIP-seq libraries for granular cells at the three stages but failed. It has been reported that PEI rewiring is often accompanied by changes in enhancer activity. Therefore, in the revised manuscript, we constructed six ChIP-seq libraries using an antibody against H3K27ac for GCs at the SWF, F1, and POF stages, with two biological duplications per stage, and annotated the genes contacting poised enhancers (PEs, 30.47%), regular enhancers (REs, 47.94%) and super-enhancers (SEs, 21.59%). These results accurately revealed the PEI regulatory networks in GCs during follicular development (Main text page 30, lines 497-507).

Comment 1-4：
"The spatial organization of compartments constructed by miniMDS 35 showed that PC1 value and gene expression are negatively associated with the distance from the center of the nucleus, i.e., the nuclear radius" -it is not clear how the nuclear radius was measured in each cell type.

Response: 1-4:
Thank you for this suggestion. We regret that the miniMDS 35 method was not described in detail in the original manuscript version. In the revised manuscript, we added this part to the Methods section:

Chromatin 3D modeling and
The 3D chromosome conformations were inferred for each Hi-C map based on the normalized intra-(at 100 kb resolution) and inter-chromosomal (at 1-Mb resolution) contact maps using an approximation of multidimensional scaling (MDS) method as implemented in miniMDS program (Rieber and Mahony, 2017). After modeling, through Euclidean distance to measure the relative distance of each chromosome (100 kb resolution) to nucleolus (start point).

Comment 1-5:
The Introduction could be improved and better structured. In the current version Introduction ends with the too long section summarizing main results. By my opinion, major goal of the study is also not clearly stated in the Introduction.

Response: 1-5:
Thank you for this comment. We have simplified the summary (as shown in the following) and summarized the main goals in the Introduction section (Lines 62-72) in the revised manuscript.

Introduction
The domestic chicken (Gallus gallus domesticus), which includes broiler (meatproducing) and layer (egg-producing) chickens, is of enormous agricultural significance and represents a classic model to study folliculogenesis (Bahr et al,. ). In this study, we investigate the transcriptomic dynamics of GCs in ovarian follicles across ten key developmental stages and generate high-resolution chromatin contact maps for GCs across three major developmental stages using in situ high-throughput chromatin conformation capture (Hi-C) sequencing. These experimental settings allowed us to conduct an integrated analysis of chromatin structure and transcriptomic characterization of chicken GCs associated with various physiological functions during folliculogenesis.

Comment 1-6:
Please compare the RNAseq data obtained for chicken ovarian follicles with previously published studies describing chicken follicle transcriptome.

Response: 1-6:
Thanks for your comment. We have discussed the results in previously published RNA-seq studies focusing on chicken ovarian follicles in the Discussion section.

Discussion
Poultry breeders have always sought chickens with high egg production, and this trait depends on efficient ovarian development and ovulation. Additional insight into the gene transcription process during folliculogenesis will help to better understand the reproductive physiology of hens and eventually improve their laying performance. Here, we performed bulk RNA-seq to systematically investigate the gene expression profiles of GCs in ovarian follicles across the whole development process, including ten stages. Substantial transcriptomic dynamics showed distinct gene expression patterns corresponding to specific stages of ovarian follicle development. We found that the SWF, F1, and POF stages, which represent the prehierarchical, preovulatory, and postovulatory phases, respectively, had the largest transcriptome differences between each other among all stages of follicle development. GO terms of stage-specific signature genes of GCs were identified by scRNA-seq and emphasized the significant functional differences in the SWF, F1, and POF stages. These changes in gene expression, such as AMH, reflect the corresponding functional characteristics of follicle development at different physiological stages, which was consistent with other RNA-seq studies in chicken follicle development (Zhu G et al., 2015;Zhu G et al., 2019).

Comment 1-7:
Please compare HiC data and identified TAD borders in chicken granulosa cells with other chicken cell types.

Response 1-7:
Thank you for your comment. To investigate the conservation of TAD in different cells, we downloaded chicken fibroblast and erithrocyte Hi-C data (including immature and mature erithrocytes) and identified TADs in these cells using methods similar to GCs. The comparison with TAD demonstrated that GCs and fibroblasts had TAD structures (mean Spearman's r of Direction Index = 0.87) that were more similar than GCs and erithrocyte cells (r < 0.32). Moreover, we identified 1996, 1361, and 1326 TAD boundaries in the chicken fibroblasts cells, immature and mature erithrocytes, of which, 82.01%, 55.84%, and 56.71% TAD boundaries are conserved in the GCs, respectively ( Figure S7 d-e).

Comment 1-8:
By my opinion, the phrase «dramatic remodeling of genome structures» is too strong for the data obtained.

Response 1-8:
Thank you for your comment. We have revised the sentence to read "Notably, we found that higher-order chromatin structures, including compartmentalization, TAD boundaries, and intra-TAD interactions, were dynamic during the stage transformation of GCs and were associated with gene expression changes" (Main text page 37, lines 594-597).

Comment 1-9:
Within the Abstract: "we integrated RNA-seq and Hi-C analyses of chicken follicular granulosa cells of 10 developmental stages." please correct the sentence so that it would be clear that HiC analyses was performed in 3 developmental stages.

Response 1-9:
Thank you for this correction. We have revised this sentence as suggested: "We investigate the transcriptomic dynamics of chicken GCs over ten follicular stages and assess the chromatin architecture dynamics and how it influences gene expression in GCs at three key stages: the prehierarchical small white follicles (SWF), the first largest preovulatory follicles (F1), and the postovulatory follicles (POF)." (Main text page 1, lines 9-13).

Response 1-11:
Thank you for this comment; we have modified the display of the Figure 4g.

Reviewer #2：
Comment 2-1: Response 2-1: We are grateful for the reviewer's help in reviewing our manuscript and for the positive feedback.

Comment 2-2:
One of the main concerns with this report is the analysis of the follicles as a whole unit, without the differentiation between the somatic and the germ compartments. The authors did not mention any limitations of the study. One limitation I could think about is the lack of the germ cell analysis.

Response 2-2:
We are grateful for your comments on our manuscript. In this study, we isolated the granulosa cells (GCs) in the follicle, which did not contain germ cells. These GCs were subsequently analyzed. We agree with the reviewer that our study is limited by a lack of germ cell analysis, which is addressed in the revised manuscript.
It is worth noting that we primarily assess GCs in this study. Advances in sequencing technology at single-cell resolution mean that further research is needed to explore the 3D genome architectural dynamics and how it influences gene expression in oocytes throughout the reproductive cycle, as well as analyze various epigenetic data to better understand regulatory mechanisms in follicle development.

Comment 2-3:
Abstract -please avoid using acronyms, such as SWF, F1, POF and TAD, in the abstract.

Response 2-3:
As suggested, we have deleted the unnecessary acronyms in both the abstract and the main text.

Response 3-1:
We appreciate the reviewer's comments and helpful suggestions for our manuscript.

Response 3-2:
We appreciate the reviewer's valuable comments. In our study, we generated an ultra-deep Hi-C contact map at 5-kb resolution (~95.13% of 5 kb bins had at least 1,000 intrachromosomal contacts) by merging the Hi-C data of two replicates. This enables us to identify PEIs at 5-kb resolution and investigate gene expression regulation (Page 29, Line 483). The deep Hi-C data can identify compartments and TADs at a 20-kb resolution, without merging contacts of replicates.
The correlation between chromatin conformation (form) and gene expression

Comment 3-3:
The authors should add a two-dimensional scatter plot to visualise the global relationships among RNA-seq samples.

Response3-3:
As suggested, we have added a scatter plot ( Figure S1b) to the revised manuscript to display the relationship between the samples, which were obtained from the t-Distributed Stochastic Neighbor Embedding (t-SNE) analysis using bulk RNA-seq data.

Response 3-4:
Thank you for your helpful comment. As suggested, the revised manuscript displays functional enrichment results in the following bubble plots: Figure 1c, Figure S1d, Figure S2f, Figure S9c, and Figure 6d in the revised version.

Comment 3-5:
Given the study's descriptive nature and lack of resolution in chromatin data, the authors should define specific active putative cis-regulatory elements by adding chromatin accessibility (ATAC-seq) characterisation at selected time points (same as Hi-C analysis). Integration of these data would greatly help understand the intricacies of the gene regulatory programme underlying chicken ovarian follicle development.

Response 3-5:
Thank you for your helpful comment. We have added the ATAC-seq analysis and included high-quality chromatin accessibility information in the revised manuscript. We also performed an ATAC-seq assay to measure differences in local accessibility during folliculogenesis. As expected, we found that the A compartments were enriched by more ATAC peaks than by their B compartments, making them more accessible. We observed that stage-specific peaks in GCs at the SWF and F1 stages are enriched in motifs corresponding to the transcription factors (TFs) in the GATA family, which are essential for the development, differentiation, and homeostasis processes. This suggests that the differentiation of GCs is highly active during these two stages (Aronson et al., 2014;Bertero et al., 2005). In contrast, POF-specific peaks in GCs are  Oncol Rep 40, 1907-1916, doi:10.3892/or.2018.6596 (2018.

Comment 3-6:
It is surprising that using pairwise differential expression analysis, the authors found "no difference in gene expression between contiguous stages of development" (line 139-140), especially between SWF/LWF, LYF/F5 and F1/POF. Increasing the sequencing depth would likely for RNA-seq libraries as well allow identification of the differentially expressed genes between these stages (see comment above).

Response 3-6:
Thank you for this helpful comment. We increased the depth of sequencing (~13.54 Gb per library) and the number of samples (six replicates per time point) (Table S1) in the revised manuscript. As expected, 39~5,580 differentially expressed genes were identified between contiguous stages. We have updated this section in the revised version (page 4 lines 96-102).
Investigating enhancer activity helps to reveal the PEI regulatory network during granulosa cell development (Main text page 30, lines 496-506).
While we did not generate follicular atlas single-cell data, we constructed singlecell libraries for granulosa cells (GCs) at the SWF, F1, and POF stages.

Fig 1b legend. It is not clear what the lines denote. It is a very informative figure
validating the temporal expression of representative genes from each cluster.
First, the clusters can be ordered by cluster 2, 1, 4 and 3 to visualise the transition more clearly. Second, the vital genes of each cluster can be noted on the side of this heatmap.

Response 3-8:
Thanks for your comment. In Figure 1b of the revised manuscript, the red lines represent mean expression levels, and the blue lines represent each gene expression in a relative cluster during development. As suggested, we have rearranged the cluster order to more clearly visualize the transition (Figure 1b).
Additionally, the genes in each cluster have been listed in Supplementary Table   S3.

Fig 2e. Authors should state the ratio of medians in addition to the p-value, as
the low p-value can be mainly caused by a large number of sample sizes.

Response 3-9:
Thank you for your comment. We included the gene number and median values in Figure. S4d and other figures throughout the revised manuscript.

Comment 3-10:
The details of multiple testing correction and other statistical methods should be clearly noted in the text, methods and figure legends.

Response 3-10:
Thank you for your comment. As suggested, we included information about the statistical test in the figure legends and the Materials and Methods section (page 49, lines 944-947).

Statistical analyses
All statistical analyses were performed by Student's t-test or Wilcoxon rank-sum test using R. *, **, and *** in the figures were represent P < 0.05, P < 0.01, and P < 0.001, respectively.

Response 3-11:
We have carefully checked and corrected typographical and grammatical errors throughout the manuscript.

Reviewers' Comments:
Reviewer #2: Remarks to the Author: The research "Dynamic transcriptome and chromatin architecture in granulosa cells during chicken folliculogenesis" reported is an important contribution to the field. Most of questions and concerns from my previous review have been addressed. I recommend to accept the manuscript with some clarifications, which are listed below.
Response 1-2: It is not clear how exactly new analysis solves the noise problem. Imagine TADs which have distance normalized bona fide (i.e. with zero noise) contact frequency e_tads, and intra-TADs having bona fide distance-normalized contact frequency e_intra. Now if we consider samples with different noise levels e_noise this will results in observed interactions (e_tads + e_noise) and (e_intra + e_noise) Could you convince the readers that changing the parameter e_noise over fixed parameters e_tads and e_intra will not affect the results of computations performed using TAD intactness of VNE math? It seems that at least TAD intactness, which is essentially ratio (e_tads + e_noise) / (e_intra + e_noise) will depend on e_noise value, and I assume that VNE will behave similarly.
Response 1-3: In fact, CTCF peaks are associated with enhancer and H3K27ac signal as well, thus H3K27ac annotation can not exclude the possibility of CTCF-mediated interactions. If annotating CTCFmediated interactions is not possible, I would recommend to highlight that some PEIs might be actually CTCF-mediated loops in results and discussion sections.
Reviewer #3: Remarks to the Author: The addition of Figure 1b-c and the discussion of different clusters significantly improved the clarity of the manuscript. The authors addressed the reviewer's concerns regarding the impurity and unintended inclusion of germ cells in the analysis. Below are a few minor typos and grammar mistakes that need to be addressed before publication. Abstract: 1. There is a typo in Line 13: "Our results provide demonstrate" -please delete the word "provide". 2. Line 15-16: Please consider replacing the word "including" with "for": "providing ample evidence for compartmentalization …" 3. Last sentence in the abstract: the authors state that their results "lay the ground work for indepth functional characterization", but it is not clear functional characterization of what?
Introduction: 1. Line 27-28: I think the more accurate statement would be: "…consists of follicles at several different developmental stages", not several follicles. 2. Lines 53-54: Please, check the grammar of the statement: "… where there is a focus on the activation of primordial follicles or growth follicle selection". What do you mean by "growth follicle selection" 3. Line 57: Typo in the word "interactions" Results: Lines 278-280: The paragraph ends with a sentence "These results support 279 the physiological course of chicken folliculogenesis." and the next paragraph starts with a similar sentence "These results provide evidence…." I am confused to which results authors refer to in the new paragraph starting on Line 280. The authors have substantially revised the manuscript entitled "Dynamic transcriptome and chromatin architecture in granulosa cells during chicken folliculogenesis" according to the reviewer's (Reviewer #3) earlier comments. The substantial revision made by authors include: added a scatter plot to visualize the global relationships among RNA-seq samples; added several bubble plots to better display both gene ontology (GO) term enrichment and statistical significance values in the exact representation; added the ATAC-seq analysis and included high-quality chromatin accessibility information; increased the depth of sequencing and the number of samples to identify DEGs between contiguous stages; generated ChIP-seq libraries using H3K27ac antibodies and investigated the enhancers (such as poised enhancers, regular enhancers, and super-enhancers) activity to reveal the promoter-enhancer interaction (PEI) regulatory network during granulosa cell development; added scRNA-seq data of granulosa cells isolated at three representative stages (SWF, F1, and POF); and improved the figure and supplementary table showing temporal expression of genes from four cluster during folliculogenesis. For some comments, particularly in comment 3-2, the authors gave sensible answers and reflected that in the revised discussion.
Overall, the revised manuscript and authors point-by-point response to Reviewer #3's earlier comments are satisfactory.
However, several additional minor issues should be corrected in the revised manuscript. Minor comments Please correct the phrases at the following sentences: Page 1, line 13. "Our results provide demonstrate the"; Page 2, line 54. "growth follicle selection" Page 4, line 98. "of the three prehierarchical stages". Is this three or four? Page 5, line 129. "post-ovulatory POF stage". Write this stage correctly (as POF stage) here and several other places.
Supplementary Fig. 1b and 1e was cited again when describing the results of scRNA-seq (line 129). Citation here is not necessary and could leads to misunderstand the results from different techniques.
"chicken granule cells" should be "chicken granulosa cells" at line 176 and many other places.
Page 5, lines 132-136. It is not clear why the chicken heart cells were used as comparative controls for scRNA-seq of GCs. Please write few lines of reason. Also, add 2-3 marker expression to confirm the identity of heart cells in suppl. fig. 2c.
Similar to above comment, explanation is needed for why the chicken fibroblasts cells and erithrocytes were used for the comparison of TADs in GCs. Correct the typo as "erythrocytes" here and several other places. Add the definition of CEF, CIME, and CME.
Please cite the figure panels in a serial manner. Supplementary Fig. 11e appears at line 499, but the panels a-d appears at line 532.

Detailed responses to reviewers
All comments provided by reviewers are given in gray italics, and our responses are in black. All revisions in the manuscript are marked in red.

Reviewer #1
Comment will depend on e_noise value, and I assume that VNE will behave similarly.

Response 1-1
We appreciate this thoughtful comment. To investigate the robustness of the TAD intactness measurement, we generated the simulated contact matrix by randomly adding signal noises to different proportions (from 5% to 50%) of the contact matrix at different noise signal levels (from 10% to 100% of the average contact frequency at each distance in the original contact matrix).
As shown in Figure S8 (newly added, Main text pages 26-27, lines 421-443), we found the TAD intactness was highly robust across variable noise levels.
Even at noise levels as high as 60% for 50% of the contact matrix, the TAD intactness between the original and simulated contact matrix is approximately equal (fold change of intactness = 1, P = 0.65, Wilcoxon rank-sum test).

Comment 1-2
In fact, CTCF peaks are associated with enhancer and H3K27ac signal as well, thus H3K27ac annotation can not exclude the possibility of CTCF-mediated interactions. If annotating CTCF-mediated interactions is not possible, I would recommend to highlight that some PEIs might be actually CTCF-mediated loops in results and discussion sections.

Response 1-2
Thank you for your constructive comments regarding the annotation of PEIs mediated by CTCF loops. Per your suggestion, we used the consensus CTCF motif of vertebrates to in silico predict the CTCF motif loci and their orientations in the chicken genome. As expected, these CTCFs were enriched in TAD boundaries and enhancer regions (newly added Figure S12k).
We also identified loops in SWF, F1, and POF at 5 kb resolution in the genomic distance range of 20 kb to 2 Mb using the Fit-Hi-C Python package.
We next applied a hard cutoff to obtain the top 15,000 loops by ranking the strength of the loops. Finally, we obtained 6,632, 7,050, and 8,023 convergent CTCF-CTCF loops in SWF, F1, and POF, respectively. The genes in these loops have 4,067, 4,926 and 2,582 PEIs respectively. As expected, the higher proportion of these PEI events (76.52%, 78.56%, and 80.38% for SWF, F1, and POF respectively) have occurred within loops (newly added Figure   S12l).
We have added a Supplementary table for these annotated convergent CTCF-CTCF loops and relative PEI information (Table S9), as well as an annotation of the CTCF loop in Fig. 6e (Main text page 33, lines 566-568).

Identification of convergent CTCF-CTCF loops
We used the FIMO software (v5.1.1) to identify CTCF motif loci and their orientations in the GRCg6a version of the chicken genome based on consensus CTCF motif from the JASPAR CORE 2016 vertebrate database. As expected, these CTCFs were enriched in TAD boundaries and enhancer regions ( Figure S12k). We also separately identified loops in SWF, F1, and POF at 5 kb resolution in the genomic distance range of 20 kb to 2 Mb using the Fit-Hi-C Python package (v2.0.7) (q value < 0.05). To further obtain the highly confident loops, we applied a hard cutoff to obtain the top 15,000 loops by ranking the loop strengths. We finally obtained 6,632, 7,050, and 8,023 convergent CTCF-CTCF loops in SWF, F1, and POF, respectively. The genes in these loops have 4,067, 4,926 and 2,582 PEIs respectively.
Newly added descriptions in Results: In addition, the convergent CTCF-CTCF loop associated with FDX1 was only identified in the F1 stage, which is consistent with the most PEIs detected in this stage (Fig. 6e). This result indicates that CTCF sites are preferentially located near promoters and/or enhancers to constrain the PEIs (Supplementary Fig. 12k-l, Supplementary Table 9).

Comment 2-1
The addition of Figure 1b