Significant heterosis for egg production in chicken reciprocal crosses
We constructed the reciprocal cross using two genetically distant chicken breeds (White Leghorn and Beijing You chicken) to decipher the molecular mechanism of egg production heterosis. The crossbreds (WY and YW) showed significant heterosis throughout the laying periods (Fig. 1a), which was divided into pre-peak laying stage from the age at first egg to 28 wks (s1), peak laying stage from 29 to 36 wks (s2), persistent laying stage from 37 to 72 wks (further split into mid laying stage from 37 to 56 wks (s3) and late laying stage from 57 to 72 wks (s4)) and prolonged laying stage from 73 to 100 wks (s5) based on the physiological status and egg production rate of hens. Compared to mid-parent value (MPV), the crossbreds showed 2.88 days (WY) and 2.23 days (YW) earlier for the age of first egg, respectively. During s1, the average heterosis for egg production was 7.44% and 7.36% for WY and YW, respectively. During persistent laying period (s2 ~ s4), the heterosis decreased in both crossbreds with the average heterosis 4.02% and 2.67%, respectively. During s5, the heterosis increased to 10.24% and 8.38% for WY and YW, respectively. Overall, the crossbreds laid 33.28 and 27.24 more eggs than the MPV.
Subsequently, the five representative laying stages, including AFE, 35, 46, 72 and 100 weeks of age were selected to harvest the ovary tissues for transcriptomic analyses. The RNA sequencing on the 120 ovarian tissues for purebreds and crossbreds (six biological replicates for each stage and each group) generated 1.31 trillion data. As shown in the Fig S1, these samples were evaluated for quality with average Q30 of 94.49%. By mapping to the chicken reference genome (GRCg7b), we obtained average mapping rate of 94.54%. Principal component analysis (PCA) showed the transcriptome profile of s5 was separated from the other 4 stages, which were interacted with each other, suggesting the relatively large transcriptomic changes in the prolonged laying period (Fig. 1b). Except for the s2, the WW and YY could be separated in each stage, while the crossbreds were not clearly clustered (Fig S2). Furthermore, the genome-wide gene expression in both purebreds and the crossbreds remained closely correlated with a value lager than 0.56 supported by the pairwise Pearson correlation coefficients (Fig. 1c), indicating highly conserved transcriptional features in the purebreds and crossbreds during the whole laying periods.
Gene regulatory networks correlate with the egg production
Using the high-quality transcriptomic data sequenced in the five laying stages, we conducted the weighted gene co-expression network analysis (WGCNA) and correlated the phenotypes with the co-expressed gene modules for both purebreds and crossbreds to identified the gene regulatory networks related to the egg production. We selected genes with MAD ranked top 50% to construct a signed, scale-free network for each genetic group. The appropriate soft-thresholding power (β) for network construction were estimated as 24, 16, 16 and 16 for WW, WY, YW and YY, respectively (Fig S3). We obtained 11,646 genes harbored in 9 modules, 11,654 genes harbored in 15 modules, 11,717 genes harbored in 16 modules and 11,718 genes harbored in 15 modules for WW, YY, WY and YW, respectively (Fig. 2a). The correlations between genes and investigated traits, including egg number (EN), oviduct length (OL) and body weight (BW) was diverse within four genetic groups. For example, the correlations in WW were largely different from that in YY (Fig. 2a), suggesting the breed difference in terms of global gene expressions and their interactions with traits. We further calculated the correlations between gene modules and traits, revealing 4 (2 positively correlated modules and 2 negatively correlated modules) and 8 modules (4 positively correlated modules and 4 negatively correlated modules) were significantly correlated with EN for WW and YY, respectively. Among these modules, the green module (MEgreen) of WW was also correlated with BW, while the MEyellow, MEpink and MEbrown were only related to EN. The red and green modules of YY were also correlated with BW, and the turquoise and blue modules of YY were related to all 3 traits (Fig. 2b). In the two crossbreds, we detected 7 modules each associated with EN (Fig. 2c). Regarding the WY, two of four positively related modules (MEred and MEgreen) were also related to BW, and the blue module was related to the three traits. For the YW, two of four positively related modules (MEmagenta and MEgreenyellow) were also correlated with BW, and all three negatively related modules (MEred, MEgreen and MEblack) were also correlated with BW.
We then tested for gene enrichment of EN-related modules against GO and KEGG pathways databases to detect the biological process (BP) and pathways involved for the four genetic groups. In total, we found that genes harbored in the EN-related modules of WW, YY, WY and YW were significantly enriched in 38, 64, 33 and 15 BPs, respectively (Table S1a). After overlapping, a total of 21 BPs, such as cell adhesion, oxidative phosphorylation, peptide metabolic process and cell migration were shared by WW and YY, while other BPs accounting for 67% of total BPs were breed-specific (Fig S4). In terms of KEGG, we identified 15, 25, 22 and 16 KEGG pathways were significantly enriched by EN-related module genes for WW, YY, WY and YW, respectively (Table S1b). By filtering pathways with gene ratio less than 5%, we found that the ribosome and oxidative phosphorylation were shared by four genetic group. Multiple pathways, such as focal adhesion, ECM − receptor interaction, vascular smooth muscle contraction and regulation of actin cytoskeleton were shared by two purebreds. The TGF − beta signaling pathway and Wnt signaling pathway were specifically identified for WY and YW, respectively. Notably, the apelin signaling pathway, cell cycle and spliceosome were enriched by module genes in the two crossbreds (Fig. 2d). These results suggested that divergent biological processes and pathways were involved in the egg production among purebreds and crossbreds, and the common networks among genetic groups were potential regulators for the formation of egg production heterosis.
Highly preserved gene regulatory networks among the crossbreds and purebreds
Given the close correlations of gene expressions and common KEGG pathways/GO terms enriched by WGCNA module genes, we assessed the conservation of gene co-expression networks between the crossbreds and the purebreds by conducting 100 times of permutation tests. By projecting the WW transcriptome dataset onto the network of YY, WY and YW, all modules were strongly preserved woth Zsummary > 10 (Fig S5a). When we took modules identified in YY as a reference network, we found that all modules in the crossbreds were at least weakly (Zsummary > 2) preserved (Fig S5b). Moreover, highly preserved gene regulatory networks were found between the two crossbreds by setting WY as a reference network (Fig S5c). The conserved gene networks were corresponded to the highly correlated transcriptional architecture underlying the transcriptome profile of ovary in the purebreds and crossbreds (Fig. 1c). Subsequently, we investigated transcriptional differences in the core gene regulatory networks, which were promisingly related to the significant heterosis for egg production in the crossbreds.
Consensus gene regulatory networks and their expression patterns between the crossbreds and purebreds
The heterosis difference (Fig. 1a) and different EN-related WGCNA modules (Fig. 2) between two crossbreds suggested that divergent regulatory networks contributed to the egg production heterosis. Therefore, we analysis the patterns of conserved gene regulatory networks for WY and YW, respectively. Based on the gene modules previously detected, we firstly constructed a conserved gene co-expression network underlying 5 laying stages among WW, YY and WY. As showing in the Fig. 3a, a total of 10,666 genes assigned to the conserved network were clustered in 14 consensus modules (CMs) consisting of 100 ~ 3,100 co-expressed genes (Fig. 3a and Table S2a), which were regulated consistently in different genetic groups. We then created CM eigengenes network and clustered CM eigengenes for each genetic group. The dendrogram showed that the CM clusters of WY were different from that of purebreds (Fig S6a). The relationships between CMs and modules previously detected showed that the modules in each genetic group (genotype-specific modules) were split into multiple CMs. For example, the EN-related green module was associated with blue, red and tan CMs in WW, and associated with brown, red and tan CMs in YY. The EN-related black module was associated with yellow, black and pink CMs in WY (Fig. 3b). We further analyzed inter-module correlations among the 13 CMs (excluding grey modules) in each genetic group by generating an eigengene network. The pair-wised correlation coefficients among CM eigengenes indicated that networks of WY were partially similar to that of two purebreds in some modules (Fig. 3c). Furthermore, we compared the expression pattern of eigengenes in the CMs, and observed differing eigengene expression patterns between crossbred and parents across 5 laying periods (Fig. 3d). For example, the blue CM eigengene expression levels increased over s1 ~ s4 and then decreased in s5 across three genetic groups, while the turquoise CM eigengene expression levels decreased over s1 ~ s4 and then increased in s5 across three genetic groups. KEGG pathway analysis indicated that genes of the blue CM were significantly enriched in ribosome, efferocytosis, MAPK signaling pathway, Wnt signaling pathway, cellular senescence, focal adhesion, etc (Fig. 3e). The genes of turquoise CM were enriched in the RNA degradation, cell cycle and ubiquitin mediated proteolysis and endocytosis were decreased during laying periods. The brown CM genes that enriched in the mRNA surveillance pathway, Notch signaling pathway and spliceosome were decreased in s1 ~ s3 and then increased in s4 ~ s5. The pink CM genes that enriched in the immune-related pathways, such as lysosome, phagosome, NOD − like receptor signaling pathway, cytokine − cytokine receptor interaction, apoptosis and toll − like receptor signaling pathway, were highly expressed in s1, s3 and s5 while lowly expressed in s2 and s4. Moreover, the genes of red CM that overrepresented in the ribosome and oxidative phosphorylation were highly expressed in s5 (Fig. 3e and Table S3a). Notably, the eigengene expression of CMs in the crossbreds was more specific compared to the purebreds, suggesting a genotype-specific perturbation of gene expression levels in these CMs (Fig. 3d).
Using same pipeline, we analyzed the transcriptional differences that are associated with the egg production heterosis among WW, YY and YW. We firstly identified 10,640 genes were cluster into 14 CMs (Table S2b and Fig S7a). Genotype-specific modules were also correlated with multiple CMs, for example, green module in YW was correlated with yellow, purple and pink CMs (Fig. S7b). The correlation heatmap among CMs showed that the expression networks of YW were specifically similar to that of two purebreds in some modules (Fig. S7c). We found CMs like blue decreased over s1 ~ s4 and then increased in s5 in three genetic groups (Fig. S7d). Genes harbored in the blue CM were significantly enriched in the cell cycle, DNA replication, progesterone-mediated oocyte maturation and so on (Fig. S7e and Table S3b). Based on results from the two crossbreds, we deduced dynamic gene expressions and interactions in various pathways were associated with egg production heterosis in the crossbreds compared with the purebreds.
Dominant gene expression for conserved biological pathways in the crossbreds
To identify key regulators that may participate in the dynamic regulation of heterosis in the crossbreds, we screened hub genes that were highly connected to other genes in the conserved gene co-expression networks for WY and YW, respectively. We calculated the module membership value (kME) for each gene, from which we identified 4,272, 4,422 and 4,217 hub genes (kME > 0.70, P < 10e-5) in WW, YY and WY, respectively (Fig. 4a and Table S4). Among these hub genes, a total of 2,392 genes were shared by three genetic groups, indicating that these genes have preserved functional features across three groups. Function enrichment analysis revealed that these shared hub genes were significantly (FDR < 0.05) overrepresented in GO biological processes related to cell cycle, DNA damage response, macromolecule catabolism and RNA catabolism (Table S5a). KEGG pathway analysis showed that the shared hub genes were enriched in the endocytosis, cell cycle, cellular senescence, FoxO signaling pathway, etc. (Fig. 4b and Table S6a), some of which were also enriched by genes in the EN-related module. As showing in Fig. <link rid="fig11">4</link>c, 1, <link rid="fig11">4</link> and <link rid="fig11">4</link> pathways were overlapped with enrichments of hub genes for WW, YY and WY, respectively. The four overlapped pathways for the WY were cell cycle, nucleocytoplasmic transport, mRNA surveillance pathway and Fanconi anemia pathway, among which common genes were found between nucleocytoplasmic transport and mRNA surveillance pathway and between cell cycle and Fanconi anemia pathway, respectively (Fig. 4d). By focusing on the 31 unique genes, we further profiled their expression patterns. The dominant expression patterns were found for genes overrepresented in the cell cycle and Fanconi anemia pathway (Fig. 4e). The cell cycle pathway regulates cell division and proliferation including cell growth and DNA replication, and Fanconi anemia pathway is made up of 22 proteins that work together to repair damaged DNA and assist with DNA replication and other cellular processes. The genes involved in the cell cycle were mostly expressed in high-parent patterns during s1, and then they were transformed into low-parent patterns during s3. The expressions of genes were higher in YY than that in WW (Fig. 5a), and the genes expressed in dominant patterns in s2, s4 and s5 were less than that in s1 and s3. Among these genes, ORC3 and CHEK2 were also dominantly expressed in low-parent pattern during s2 and s4, respectively. CDK7 and MAD2L1 were expressed in below-parent pattern, while MTBP and ATR were in low-parent pattern. During s5, OCR3 and ATR were expressed in below-parent pattern and high-parent pattern, respectively. CHEK2 encodes a protein that functions as a regulator of the cell cycle as well as a tumor suppressor. The protein is activated in the presence of DNA damage to prevent the entry into mitosis. MAD2L1 accumulates early in oocyte maturation and affects cell proliferation and cell cycle progression.
In the combination of WW, YY and YW, we identified 1,866 common hub genes, which were enriched in the GO BPs related to cell cycle and macromolecule catabolic process (Table S5b) and KEGG pathways including efferocytosis, cellular senescence, cell cycle and FoxO signaling pathway (Fig. S8 and Table S6b). Among these pathways, cell cycle was the only KEGG pathway overlapped with pathways enriched by EN-related module genes (Fig. 4c). Again, we profiled expression patterns for 6 common pathway genes and found that the 4 identical genes including MAD2L1, CDC16, MTBP and ATR were additively expressed or dominantly expressed in high-parent model during s3 (Fig S8c). This is probably the underlying factors leading to the heterosis difference between 2 crosses. Another 2 genes including PRKDC and E2F1 were low-parent expressed during s4 and during s3 and s5, respectively. Similar low-parent expression pattern between WY and YW for MTBP and ATR were found during s4, while the expressions were higher in WW than that in YY (Fig. 5b). E2F1 was one of transcription factors of E2F family, playing an important role in the control of the cell cycle, cell proliferation, and differentiation. PRKDC encodes the DNA-dependent protein kinase catalytic subunit (DNA-PKcs) protein, which plays an important role in nonhomologous end joining of DNA double-strand breaks and is also closely related to the establishment of central immune tolerance and the maintenance of chromosome stability. These results supported the key roles of hub genes involved in the cell cycle pathway played in egg production heterosis during 5 laying stages in the reciprocal crosses.