- Research article
- Open Access
Multi-omics analysis reveals that natural hibernation is crucial for oocyte maturation in the female Chinese alligator
BMC Genomics volume 21, Article number: 774 (2020)
Hibernation in an appropriate environment not only is important for the survival of hibernators in winter, but also is crucial for breeding in the following season for many hibernating species. However, the genetic and epigenetic mechanism underlying this process remain unclear. In the current study, we performed an integrative multi-omics analysis of gonads collected from Chinese alligators that overwintered in wild cave and artificial warmroom to explore transcriptomic and epigenomic alternations in these organs.
The data revealed that in the breeding season, female alligators were more strongly affected in terms of gene expression than males by non-hibernation because of overwintering in a warm room, especially for genes related to oocyte maturation, and this effect commenced in winter with the downregulation of STAR, which is the rate limiting factor of steroid biosynthesis. Further, miRNAs were found to play essential roles in this negative effect of overwintering in the warm room on hibernation. The upregulated miRNAs likely were responsible for the suppression of oocyte maturation in the breeding season. Finally, DNA methylome changes, especially hypomethylation, were found to play an important role in the alterations in ovarian function-related gene expression induced by non-hibernation.
Our study revealed the crucial role of hibernation quality for oocyte maturation in the Chinese alligator and the underlying genetic and epigenetic mechanisms, and highlights the importance of habitat, and especially, the overwintering site, in the conservation of not only the Chinese alligator, but also other endangered hibernators.
Hibernation is an adaptive strategy adopted by many animals to survive the cold and foodless winter . During winter, hibernators stay in their refuges and suppress their metabolic rates to substantially save energy, water, and oxygen . The underlying molecular mechanisms vary between endotherms and ectotherms [3, 4]. Increasing evidence gathered since mid-last century suggests that there is a close correlation between the overwintering patterns and reproductive cycles of many mammals [5,6,7], birds , reptiles , amphibians , and invertebrates . However, the importance of hibernation in gonadal development and reproduction remains controversial. For example, hibernation is crucial for oviposition and fertilization in the mountain yellow-legged frog (Rana muscosa) . In the boreal toad (Bufo boreas boreas), hibernation would improve the breeding success in males, but not females . However, a recent study on the Wyoming toad (Anaxyrus baxteri) revealed that hibernation is not essential for optimal sperm metrics in this species . Thus, the relation between hibernation and reproduction varies among species. Further, our understanding of the molecular mechanisms underlying the effects of hibernation on reproduction is limited.
The Chinese alligator is a critically endangered freshwater crocodilian that once was widely distributed in China. However, because of climate change, hunting practices, and habitat destruction, its natural habitat is currently restricted to the lower Yangtze River, and the majority of animals live in captivity in two nature reserves in the Zhejiang and Anhui provinces of China. Chinese alligators stop eating and go into hibernation when the air temperature decreases in late October, and they emerge in late March as the temperature rises, and mate in late May [9, 14]. Hibernation quality has a crucial impact on their health and reproduction, and alligators made to overwinter in a warm room do not hibernate. When Chinese alligators move to a warmer area or overwinter in the discomfort of an artificial environment, they become less fertile [9, 15,16,17]. Thus, the Chinese alligator can serve as an ideal model to study the role of hibernation in reproduction, and the underlying molecular mechanisms.
Several genes involved in ovarian development and breeding activity, including ESR1 (encoding estrogen receptor α) , FSH and FSHR (encoding follicle-stimulating hormone and its receptor) [19, 20], and Kiss1 and Kiss1R (encoding Kisspeptin1 and its receptor) , are upregulated in the ovaries and several other tissues during the breeding season as compared to hibernation in the Chinese alligator. These results indicate that these genes may play a role in ovarian development and breeding activity after waking from hibernation; however, current knowledge is limited to a limited number of genes. Seasonal transcriptome changes in non-gonad tissues have been investigated in the Chinese alligator [4, 22] and other hibernators [3, 23,24,25,26,27,28]. We reasoned it would be interesting to use mRNA sequencing (mRNA-seq) to characterize genome-wide regulatory networks underlying the effects of hibernation on breeding success.
DNA methylation and small RNAs respectively are pre- and post-transcriptional epigenetic regulatory mechanisms. In this study, we carried out mRNA-seq, small RNA sequencing (sRNA-seq), and bisulfte sequencing (BS-seq) analyses of gonads collected in different seasons from Chinese alligators that overwintered in different overwintering environments to explore transcriptomic and epigenomic alterations induced by lack of hibernation.
We carried out RNA-seq, sRNA-seq, and BS-seq analyses of ovary and testis samples collected in winter and summer from individual alligators that overwintered in the wild (hibernation) or in a warm room (non-hibernation) to gain insights into the molecular mechanisms underlying the effects of hibernation quality on gonadal development (Fig. 1a).
The eight stand-specific mRNA-seq libraries generated a total of 468.72 million reads, 419.72 million (89.55%) of which were uniquely mapped to the Chinese alligator reference genome . mRNA-seq data analysis was based on the mapped reads, using a Chinese alligator gene annotation database containing 27,500 protein-coding genes . sRNA-seq of the eight sRNA-seq libraries generated 100.63 million clean single-end reads, 85.63 million (85.10%) of which were 18–35 nt in length and were uniquely mapped to the Chinese alligator genome. sRNA-seq data analysis was based on the mapped reads and Chinese alligator sRNA annotation data generated in our previous study . The eight BS-seq libraries generated 2955.58 million clean read pairs (669.17 Gb clean data), with an average depth of 18.38 per strand for each sample. On average, 92.65% of cytosines were covered by at least five unique reads in each sample. The bisulfite conversion rate was high in all libraries (99.51% on average).
Transcriptional alterations in the maldeveloped ovary
To obtain an overview of the gonadal transcriptome changes in female and male alligators overwintering without hibernation, we carried out principal component analysis (PCA) of protein-coding gene and miRNA expression data from the eight gonadal samples. The results suggested that overwintering in the warm room led to ovarian maldevelopment, whereas testes were less affected. The four testis samples clustered together, and samples collected in the same season were located closer to each other in the PCA plot. The ovary sample collected in summer from the individual overwintering in the warm room (SF_OVA_R) was clearly separated from the three other ovary samples (Fig. 2a and b). Pearson correlation coefficients of mRNA-seq and sRNA-seq data for gonad samples were consistent with the PCA results (Fig. 2c and d). Comparison of gene expression levels in samples collected from animals overwintering in the warm room (WM_OVA_R, WM_TES_R, SF_OVA_R, and SF_TES_R) with those in corresponding samples of hibernating or hibernated animals (WF_OVA, WM_TES, SF_OVA, and SM_TES) revealed that there were substantially more differentially expressed genes (DEGs) in ovaries collected in summer (SF_OVA_R vs. SF_OVA) than in the other three pairs. Whereas there were only 474 and 358 DEGs in the testes collected in winter and summer and 210 DEGs in the ovaries collected in winter, there were 3412 DEGs in the ovaries collected in summer (Fig. 2e). These results suggested that hibernation is more important for ovarian than for testicular development and that gene/miRNA expression may substantially alter in the ovaries during the breeding season if the animal cannot hibernate in an appropriate environment. This is consistent with a previous report that overwintering in a man-made cave or in the warmer southern region of China led to maldevelopment of the gonads and breeding failure in the Chinese alligator . Our results also suggest that this failure is likely caused by ovarian rather than testicular disorder.
Next, we investigated KEGG pathway enrichment of DEGs that were downregulated in the SF_OVA_R sample as compared to the three other ovary samples. ‘Cell cycle’ and ‘oocyte meiosis’ were enriched in DEGs suppressed in the SF_OVA_R sample (q < 0.1) (Additional file 1: Table S1), suggesting that the breeding failure caused by overwintering in an inappropriate environment is due to suppression of oocyte meiosis and maturation. Indeed, a larger number of genes in the oocyte meiosis pathway were significantly suppressed in the SF_OVA_R sample (Wilcoxon signed-rank test, q < 0.05) (Fig. 3a and b). In particular, maturation-promoting factor (MPF), the key driver of meiotic progression, was suppressed in the SF_OVA_R sample. Genes in the Mos-MEK1-ERK2, Plk1-Cdc25, Bub1-Mad1/2-APC/Cdc20, and Camk2-Emi2 pathways, through which progesterone and fertilization induce resumption of the two meiotic division cycles and oocyte maturation, were also downregulated (Fig. 3c). These results are consistent with our finding that the ovarian follicle development is depressed in SF_OVA_R compared with that in SF_OVA (Fig. 1b). Several genes crucial for fertility (GDF9, PRLRB, ADAMTS1, RBP4, and NANOS1) [30,31,32,33,34], ovarian development (BMP2A and BMP2B) , and steroid hormone biosynthesis (STAR, CYP11A, HSD17B1, LHX9, WT1, and SF1) [36, 37] were also suppressed in the SF_OVA_R sample (Additional file 2: Figure S1). Interestingly, several male sex determination factors (including MAP 3 K1, PTCH1, SOX9, FGFR2A, and AR) [36, 37] were upregulated in this sample, suggesting that repression of the male gonadal developmental pathway was disturbed (Additional file 2: Figure S1). Together, these data indicated that gene expression in the ovaries is severely distorted in the breeding period when the Chinese alligator cannot hibernate in an appropriate environment.
Although females overwintering in the warm room did not show significant transcriptome alterations in the ovaries until the breeding period, changes in the ovaries had been triggered in winter. STAR, which encodes the rate-limiting enzyme steroidogenic acute regulatory protein, was downregulated in the ovaries of alligators overwintering in the warm room (WF_OVA_R) when compared with the levels in alligators hibernating in the wild (WF_OVA) (Additional file 2: Figure S1a and c). In contrast, FGFR2A, which regulates proliferation and Sertoli cell differentiation during testicular development, was upregulated in the WF_OVA_R sample, suggesting that the higher temperature in the warm room revoked the repression of male gonadal development during winter (Additional file 2: Figure S1a and c).
MiRNA transcription alterations in the maldeveloped ovary
To explore the roles of miRNAs in non-hibernation-induced infertility in the Chinese alligator, we constructed a heat map of the expression levels of differentially expressed miRNAs (DEmiRs) among the four ovary samples, and we classified all 96 DEmiRs into 12 groups based on their expression pattern (Fig. 4a). MiRNA expression clustering revealed a clear separation of the SF_OVA_R sample from the other three ovary samples, which was consistent with the PCA results (Fig. 2b). The expression patterns of these miRNAs suggested their potential involvement in ovarian development and function. For example, in the breeding season, miRNAs in group A2 were enhanced (e.g. miR-128 and let-7c) in the hibernated individual, but not in the individual that overwintered in warm room, suggesting that they may play crucial roles during the breeding season in the Chinese alligator and that overwintering in inappropriate environment inhibits their upregulation. MiRNAs in group C2 (e.g. miR-9a and miR-790) may be critical for ovarian development and function in both seasons as they their expression levels remained high in hibernating/hibernated animals, but decreased in individuals that spent their winter in the warm room. In contrast, miRNAs in groups B3 (e.g. miR-449a and miR-135) and D3 (e.g. miR-32 and miR-153) were upregulated during summer and winter respectively in the individuals overwintering in warmroom, suggesting their potential harmful impact on ovarian development and function.
To investigate the roles of miRNAs in gene expression regulation in non-hibernation-induced infertility in the Chinese alligator further, we predicted the targets of each miRNA and generated a regulatory network of their involvement in oocyte maturation (Fig. 4b). The results showed that most of the suppressed genes related to oocyte maturation were targeted by miRNAs (Fig. 4b) and, as expected, most of these miRNAs belonged to group B3 (Fig. 4a). For example, while miR-200b was upregulated, transcription of its target MOS was substantially supressed in the SF_OVA_R sample. Similar patterns were observed for other miRNA–gene pairs, including miR-129a–ACP12, miR-210–CDK1, and miR-449a–PLK1 (Fig. 4c). These results suggest that oocyte maturation is actively repressed in the breeding season, at least in part, via miRNAs when the Chinese alligator did not hibernate in a suitable environment. In addition, the downregulation of miRNAs targeting male sex-determining genes (such as PTCH1 and AR) likely results in the upregulation of these genes’ expression, which disturbs ovarian development in the breeding season (Fig. 4c). These results suggested that miRNAs play essential roles in ovarian maldevelopment induced by overwintering in an artificial environment.
DNA methylome changes in the maldeveloped ovary
To study DNA methylation changes in ovarian maldevelopment of the Chinese alligator, we compared the pairs of ovary samples and identified differentially methylated regions (DMRs) in the genomes. We found that an average of 26.13 Mb (1.186% of the genome) were DMRs between the SF_OVA_R sample and the three other ovary samples, and an average of 2.14 Mb (0.097% of the genome) were DMRs between each pair of the three other ovary samples (Fig. 5a), whereas only 1.51 Mb (0.069% of the genome) were DMRs among the four testis samples (Fig. 5b). These patterns were consistent with the transcriptome data (Fig. 2a, b).
We next focused on DMRs between the two ovary samples collected in summer to explore the epigenetic mechanism underlying ovarian maldevelopment. Compared with the ovary sample of the hibernated alligator (SF_OVA), there were 26.03 Mb DMRs in the ovary sample of the non-hibernated alligator collected (SF_OVA_R), including 20.48 Mb hypo-DMRs (corresponding to 5679 differently methylated genes, DMGs) and 5.54 Mb hyper-DMRs (corresponding to 2304 DMGs). The overlap of the DEGs and DMRs revealed that DMGs were likely to be DEGs (Fisher exact test, q = 1.079E-148) (Additional file 3: Fig S2a and b). In the SF_OVA_R sample, genes with hypo-DMRs in their promoters or/and gene bodies were more likely to be upregulated DEGs (Fisher’s exact test, q < 0.05) (Additional file 3: Fig. S2c and d), whereas only genes with hyper-DMRs in their promoter tended to be downregulated (Fisher’s exact test, q < 0.05) and those with hyper-DMRs in their gene body did not show significant trends (Fisher’s exact test, q = 0.203) (Additional file 3: Fig. S2e and f). In addition, the average DNA methylation levels across gene elements of upregulated DEGs were lower in SF_OVA_R than in the other samples (paired-samples T test, q < 0.05) (Fig. 5c–f). These results suggested that DNA methylome changes, especially hypomethylation, play an important role in gene regulation in ovarian maldevelopment caused by non-hibernation in winter.
To investigate the regulatory role of DNA methylation in ovarian maldevelopment further, we compiled a list of 157 genes that likely play roles in oocyte meiosis, female fertility, and sex determination, and investigated their methylation status (Additional file 4: Table S2). In total, 69 of these genes (45.586%) displayed significantly differential methylation patterns between SF_OVA_R and SF_OVA (Fig. 5g and Additional file 4: Table S2). As expected, there were many more hypo- than hyper-DMGs in summer (38 vs 12), and 19 genes contained both hyper- and hypo-DMRs (Fig. 5g and Additional file 4: Table S2). Among the nine DEGs with DMRs in their promoters, seven genes showed opposite changes in expression and methylation alterations. For example, the male sex determination genes SOX9 and PTCH1, upregulation of which likely contributed to ovarian maldevelopment, had hypo-methylated promoters in the SF_OVA_R sample. In contrast, STAR was repressed and had a hyper-DMR in its promoter. However, there were more numerous DMGs with DMRs in their gene body, and the correlations between their methylation and expression were very complex. For example, some of the genes related to oocyte meiosis and downregulated in the SF_OVA_R sample had hyper-DMRs in their gene body (CYCB), some had hypo-DMRs in their gene body (AURKA, PLK1, CDC25, YWHAG_H and BUBL1), and others had both (CPEB2 and ITPR1) or none (Fig. 5g).
Seasonal transcriptional and epigenetic alterations in the testes
As mentioned above, the overwintering environment did not have a significant influence on seasonal transcriptional and epigenetic alterations in the testes. KEGG pathway and GO term analyses revealed that KEGG pathways involved in primary metabolic processes, gene expression, and protein metabolic processes, including ‘ribosome’, ‘spliceosome’, ‘RNA polymerase’, and ‘pyrimidine metabolism’ (q < 0.05, Additional file 5: Table S3) and the GO terms ‘ribosome’, ‘translation’, ‘gene expression’, ‘primary metabolic process’, and ‘protein metabolic process’ were enriched in DEGs upregulated in the testis samples of hibernating animals (q < 0.05, Additional file 6: Table S4). These results were in contrast with gene expression patterns reported in other tissues, in which energy metabolism and functional gene expression were suppressed , suggesting that the testes of the Chinese alligator actively develop during winter in preparation for reproduction in the next summer, but this process is not effect by the quality of overwinter environment. AMH (encoding anti-Müllerian hormone), MTNR1A (encoding melatonin receptor 1A), and the steroid hormone biosynthesis genes CYP11A, HSD3B, and CYP17 were upregulated in summer (Additional file 7: Figure S3), suggesting the important roles of anti-Müllerian hormone, melatonin, and steroid hormone in male breeding.
Hibernation quality is known to have a substantial impact on reproduction in the Chinese alligator [9, 15]. In this study, we found that not only the transcriptome, but also the miRNome and DNA methylome are more strongly affected in female Chinese alligators by an inappropriate hibernation environment than in male animals. DEGs, DEmiRs, and DMGs potentially involved in this process were identified to uncover the complex molecular mechanisms underlying the influence of hibernation on the breeding success of the Chinese alligator.
While hibernation quality is known to affect animal reproduction in various species, its effects seem to vary among species and the sexes [10, 12, 13]. In the testes of the Chinese alligator, the transcriptome and epigenome showed season-biased changes (e.g. in hormone biosynthetic genes), but were not obviously affected by the overwintering environment (Fig. 2 and Fig. 5b). However, in the ovary, mRNA and miRNA expression patterns as well as DNA methylation patterns were significantly altered in the breeding season following a non-hibernation winter in the warm room (Fig. 2 and Fig. 5a). These results suggest that the female alligator is more sensitive to the hibernation environment. In line herewith, it is believed that the low temperatures experienced during hibernation are an important environmental factor for oocyte maturation in frogs and other hibernating species . In the Chinese toad (Bufo Bufo gargarizans), low temperatures during hibernation contribute to oocyte maturation competence via the accumulation of CDC2, an MPF component, in the oocytes, and the CDC2 level was significantly reduced in toads disturbed from hibernation as a result of high temperatures . Our results are consistent with these findings in Anura, and in the Chinese alligator ovary, in addition to MPF components, the entire oocyte maturation pathway was suppressed when an appropriate overwinter environment was lacking.
The fact that miRNAs play key roles in gonadal development is supported by knockout/knockdown studies and expression studies in mammals [40, 41]. Knockout or knockdown of miRNA machinery genes (Dicer, Ago2, and DGCR8) in mice led to a global suppression of miRNA functions, suppression of oocyte maturation, and abnormal preimplantation development [42, 43]. In our study, miRNA expression profiling identified 96 miRNAs that were differentially expressed between the four ovary samples and that likely played roles in ovarian maldevelopment in non-hibernation (Fig. 4a). Some of these miRNAs have been reported to be related to ovarian function. For example, miR-126, miR-15b, miR-143, and miR-199a, which are abundant in the mammalian ovary and are predicted to target ovarian function-related genes, were supressed in the breeding season after overwintering in the warm room. In addition, numerous miRNAs involved in oestradiol inhibition  were upregulated during winter (groups D1–D3: miR-148a, miR-101, miR-146, miR-7b, miR-34a, miR-19, and miR-32) or summer (groups B1–B3: let-7 g, miR-144, miR-103b, miR-92a, miR-10a, miR-135, miR-129a, miR-31) in female alligators that overwintered in the warm room, suggesting that oestrogen hormone biosynthesis in the ovaries is suppressed in female alligators that missed hibernation. Consistent herewith, STAR, which encodes a rate-limiting enzyme in the steroid biosynthesis pathway, was significantly downregulated in winter in the non-hibernating female alligator.
DNA methylation is crucial in ovarian development in various species, including the honeybee (Apis mellifera) [45, 46], Japanese flounder (Paralichthys olivaceus) [47, 48], zebrafish (Danio rerio)  and sheep (Ovis aries) . Further, it participates in diseases that are of ovarian origin and directly cause infertility . We found numerous DMRs, especially hypomethylated regions, in DNA collected in summer from female animals that missed hibernation, and we identified several DMGs related to ovarian function. Moreover, as the DMGs were enriched in TF genes (Fisher’s exact test, p = 6.386E-021), and given the trans-regulatory role of methylation in hibernation adaption of non-gonadal tissue , we reason that other DEGs in the SF_OVA_R sample were trans-regulated by the differentially methylated TFs as shown in the transcription network in Additional file 8: Figure S4.
In face of the threats from climate change, hunting, and habitat destruction, the wild Chinese alligator population has fallen sharply. Therefore, it is imperative to establish a captive population and implement a captive-breeding program for population rejuvenation. However, as the quality of the hibernation environment has a great impact on oocyte maturation in this alligator, the captive environment, and especially, the potential hibernation sites, should be paid more attention to and should more closely reflect overwintering sites in their wild habitat. In addition, human disturbance should be avoided in the hibernation season.
Overall, our data provide genetic and epigenetic insights into the crucial significance of natural hibernation for reproduction of the Chinese alligator, especially for female oocyte maturation, and are expected to facilitate the development of scientific programs for successful conservation of not only the Chinese alligator, but also other endangered hibernators.
Sample sources and DNA and RNA extraction
Chinese alligator genital tissues were provided by the Changxing Yinjiabian Chinese Alligator Nature Reserve (CCANR) (Fig. 1a). Winter testis (WM_TES) and ovary (WF_OVA) samples were collected from two adult hibernating Chinese alligators (one male and one female) in January 2015 (the coldest time in Changxing), and summer samples (SM_TES and SF_OVA) were collected from another two adult individuals (one male and one female) during the breeding season, in May 2015. According to the guidelines of the Animal Ethics Committee of Zhejiang University, these animals were euthanized under deep anaesthesia by injecting ketamine (5–10 mg/kg) into the tail muscle and then followed by blood-letting. We collected not only gonads but also several other tissues used in another study. Another four adult individuals (two males and two females) were housed in a glass warm room at the CCANR in October 2014, where the alligators keep active during winter. Testis (WM_TES_R) and ovary (WF_OVA_R) samples were surgically obtained from one male and one female under deep anaesthesia by injecting ketamine (5–10 mg/kg) in January 2015. The other two individuals were released into the wild in late March 2015, and genital tissues (SM_TES_R and SF_OVA_R) of these animals were collected in the breeding season. The warmroom-overwintering alligators were released to the Changxing Yinjiabian Chinese Alligator Nature Reserve after their wounds healed well. The gonads are located in the posterior part of the abdominal cavity, in front of the kidneys. The testes are about 14*3*2 cm3 in size. The mature ovary is about 20*6*3 cm3, with many golden yellow follicles on it; while the undeveloped or maldeveloped ovary is about 10*3*2 cm3, on which the follicles are much smaller (Fig. 1b).
All samples were immediately stored in liquid nitrogen until use. All the experimental procedures with Chinese alligator in this study had been given prior approval by the Animal Ethics Committee of Zhejiang University (ZJU2015–154-13) and the State Forestry Administration of China [Forest Conservation Permission Document (2014) 1545].
The gDNA used for BS-seq and RNA used for mRNA-seq were extracted from the samples using an AllPrep DNA/RNA Mini Kit (Qiagen, Hilden, Germany), according to the manufacturer’s instructions. Total RNA for sRNA-seq was extracted using a TRIzol RNA isolation kit (Invitrogen, Waltham, MA) according to the manufacturer’s instructions.
Strand-specific cDNA library construction and sequencing
Strand-specific cDNA libraries were generated from 3 mg RNA per sample, using a NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (New England Biolabs, Ipswich, MA) according to the manufacturer’s instructions. The index-coded samples were clustered using a TruSeq PE Cluster Kit v3-cBot-HS on a cBot Cluster Generation System (Illumina, San Diego, CA). The libraries were sequenced on an Illumina HiSeq 2500 platform, generating 125-bp paired-end reads. Reads containing adapter or poly-N sequences, and low-quality reads were removed using in-house Perl scripts, and all subsequent analyses were based on the clean reads.
sRNA library construction and sequencing
sRNA libraries were prepared from 3 mg RNA per sample, using a NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB) according to the manufacturer’s instructions. The index-coded samples were clustered using a TruSeq SR Cluster Kit v3-cBot-HS on a cBot Cluster Generation System. After cluster generation, the sRNA libraries were sequenced on an Illumina HiSeq 2500 platform, generating 50-bp single-end reads. Reads containing 5′ adapter or poly-N sequences, reads without 3′ adapter or the insert tag, and low-quality reads were removed using in-house Perl scripts, and all subsequent analyses were based on the clean reads.
BS library construction and sequencing
The gDNA was spiked with non-methylated λ DNA fragments for bisulfite conversion quality control purposes. For each sample, 6 mg gDNA and 30 ng λ DNA were mixed and fragmented into 200–300 bp using sonication. After end repair and acetylation, barcodes with methylated cytosines were added to the fragmented DNA. DNA bisulfite conversion was carried out twice using an EZ DNA methylation-Gold™ Kit (Zymo Research, Irvine, CA). The DNA fragments were then amplified by PCR using KAPA Hifi HotStart Uracil+ ReadyMix (Kapa Biosystems, Wilmington, MA). The index-coded samples were clustered using a TruSeq PE Cluster Kit v3-cBot-HS on a cBot Cluster Generation System (Illumina) according to the manufacturer’s instructions. The BS library was sequenced on the Illumina HiSeq 2500 platform, generating 125-bp paired-end reads. Low-quality reads, reads containing adapter or poly-N sequences, and reads shorter than 36 nt following adapter removal were filtered out using in-house Perl scripts, and all subsequent analyses were based on the clean reads.
RNA-seq data analysis
An index of the Chinese alligator reference genome was built using Bowtie2 , and the clean reads were mapped to the reference genome using TopHat v. 2.0.12 . Reads mapped to each gene were counted using HTSeq v. 0.6.1  and differential gene expression between sample pairs was analysed using the DEGseq R package v. 1.12.0 . Genes with a false discovery rate (FDR) < 0.05 and a |log2fold-change| > 1 were considered DEGs. The number of fragments per kilobase of exon per million (FPKM) mapped fragments was calculated to estimate gene expression levels. Gene expression level heatmaps were constructed using TBtools . The DESeq R package  was used to carry out PCA of mRNA-seq data of the gonad samples and for plot construction.
sRNA-seq data analysis
Clean sRNA reads were mapped to the Chinese alligator genome reference sequence using Bowtie v. 2.2.3  without any mismatch allowed. We then annotated the reference genome sequence using the Rfam database  and RepeatMaker . Reads corresponding to rRNAs, tRNAs, snRNAs, snoRNAs, repeat sequences, exons, and introns were removed. miREvo v. 1.1  and mirdeep2  were used to predict miRNAs through exploration of the secondary structure, the Dicer cleavage site, and the minimum free energy of the reads. Predicted miRNAs were subjected to Rfam  for miRNA family analysis and identification of conserved miRNAs (found in at least one other species). In-house scripts were used to obtain miRNA counts. miRNA expression levels were normalized as the number of transcripts mapped to the miRNA per million transcripts (TPM). The DEGseq R package was used to analyse differential miRNA expression between paired samples . miRNAs with an FDR < 0.05 and a |log2fold-change| > 1 were regarded DEmiRs. MiRanda, RNAhybrid, and PITA were used to predict target genes of the miRNAs [62,63,64]. Target genes approved by at least one software package were considered targets. The DESeq R package  was use to carry out PCA of sRNA-seq data of the gonad samples and for plot construction.
BS-seq data analysis
BS-seq reads were aligned to the Chinese alligator reference genome using Bismark (v. 0.12.5) with default parameters . The reference genome was transformed into fully BS-converted versions termed ‘T genome’ (C-to-T converted) and ‘A genome’ (G-to-A converted) and then indexed using Bowtie2 . All cytosines of the BS-converted reads were transformed to thymines, and the reads were aligned to the ‘T genome’. All guanines of the BS-converted reads were transformed to adenosines and the reads were aligned to the ‘A genome’. The sequence reads that produced a unique best alignment from the two alignments (original Watson and Crick strand) were aligned back to the original reference genome to infer the methylation state of all cytosines in the sequence reads. Multiple reads mapped to the same genome regions were regarded clonal duplicates and were removed to avoid inaccuracy that might be caused by PCR amplification bias. The BS library conversion rate was estimated as the percentage of cytosines sequenced at cytosine positions in the λ-DNA reference genome.
To identify methylation sites, we modelled the sum mC of methylated counts as a binomial (Bin) random variable with methylation rate (r), as mC ~ Bin (mC + umC*r).
The methylation level of each cytosine site was quantified as the number of reads containing an mC at the site of interest divided by the total number of reads covering the cytosine site. The methylation level of a specific region was quantified as the average methylation level of all cytosine sites in this region.
DMRs between two samples were identified using swDMR . As most of the mC were in the CG context and the methylation levels in non-CG contexts were very low, we focused solely on CG sites. The sliding window was set to 1000 bp with a step length of 100 bp. To ensure statistical power, only windows with at least 10 CG sites and a coverage of at least 5 in each of the two compared samples were considered. Fisher’s exact test was used and only windows with FDR-adjusted P < 0.05 and a greater than twofold methylation level change were considered DMRs. Genes containing DMRs in their putative promoter or/and gene body regions were regarded DMGs.
GO and KEGG enrichment analyses
GO term and KEGG pathway enrichment analyses were carried out using the GOseq 2.12 package and the KOBAS package, respectively [67, 68]. GO terms and KEGG pathways with an FDR (q value) < 0.05 were regarded significantly enriched.
Availability of data and materials
The Chinese alligator reference genome is available from GenBank (assembly accession: GCA_000455745.1). The BS-Seq, mRNA-Seq and sRNA-Seq data generated in this work have been deposited in the NCBI SRA database under BioProject accession numbers PRJNA556094, PRJNA556093, and PRJNA556092, respectively.
Changxing Yinjiabian Chinese alligator nature reserve
Overwinter in warm room
Small RNA sequencing
Differentially expressed gene
Differentially expressed miRNA
Differentially methylated region
Differently methylated gene
Principal component analysis
False discovery rate
The number of fragments per kilobase of exon per million
The number of transcripts mapped to the miRNA per million transcripts
Estrogen receptor α
Follicle-stimulating hormone receptor
steroidogenic acute regulatory protein
Melatonin receptor 1A.
Yahner RH. Winter Strategies. In: Wildlife behavior and conservation. New York: Springer; 2012. p. 139–43.
Staples JF. Metabolic flexibility: hibernation, torpor, and estivation. Compr Physiol. 2016;6(2):737–71.
Faherty SL, Villanueva-Canas JL, Blanco MB, Alba MM, Yoder AD. Transcriptomics in the wild: hibernation physiology in free-ranging dwarf lemurs. Mol Ecol. 2018;27(3):709–22.
Lin JQ, Huang YY, Bian MY, Wan QH, Fang SG. A unique energy-saving strategy during hibernation revealed by multi-omics analysis in the Chinese alligator. iScience. 2020;23:101202.
Barnes BM. Relationships between hibernation and reproduction in male ground squirrels. In: Adaptations to the cold: tenth international hibernation symposium. Armidale: University of New England Press; 1996. p. 71–80.
Morrow GE, Jones SM, Nicol SC. Interaction of hibernation and male reproductive function in wild Tasmanian echidnas Tachyglossus aculeatus setosus. J Mammal. 2016;97(3):852–60.
Mendonca MT, Hopkins WA. Effects of arousal from hibernation and plasma androgen levels on mating behavior in the male big brown bat, Eptesicus fuscus. Physiol Zool. 1997;70(5):556–62.
Goodenough AE, Coker DG, Wood MJ, Rogers SL. Overwintering habitat links to summer reproductive success: intercontinental carry-over effects in a declining migratory bird revealed using stable isotope analysis. Bird Stud. 2017;64(4):433–44.
Chen BH, Hua TM, Wu XB, Wang CL. Research on Chinese Alligator. Shanghai: Shanghai Scientific and Technological Education Publishing House; 2003.
Santana FE, Swaisgood RR, Lemm JM, Fisher RN, Clark RW. Chilled frogs are hot: hibernation and reproduction of the endangered mountain yellow-legged frog Rana muscosa. Endanger Species Res. 2015;27(1):43–51.
Olsson TI. The effect of wintering sites on the survival and reproduction of Gyraulus-Acronicus (Gastropoda) in a partly frozen river. Oecologia. 1988;74(4):492–5.
Roth TL, Szymanski DC, Keyster ED. Effects of age, weight, hormones, and hibernation on breeding success in boreal toads (Bufo boreas boreas). Theriogenology. 2010;73(4):501–11.
Poo S, Hinkson KM, Stege E, Terrell KA. Sperm output and body condition are maintained independently of hibernation in an endangered temperate amphibian. Reprod Fert Develop. 2019;31(4):827.
Fang LM, Zhai T, Zhao L, Fang SG. Deep brumation features of Zhejiang Chinese alligators. Chin J Wildl. 2015;36(3):284–7.
Xia TS, Zhou KH, Zhu JL. Influence of Hibernant temperature on reproduction of Chinese alligator. Sichuan J Zool. 2006;25(2):298–399.
Xia TS, L. ZJ, M S: relationship between temperature and diseases occurrence of Chinese alligator during hibernation. Sichuan J Zool. 2006;25(2):400–2.
Zhang GL, Geng YJ, XIao JG, Yang SH. Comparison of two overwintering ways for Chinese alligator in captivity. J Econ Anim. 2003;7(2):57–9.
Zhang R, Hu Y, Wang H, Yan P, Zhou Y, Wu R, Wu X. Molecular cloning, characterization, tissue distribution and mRNA expression changes during the hibernation and reproductive periods of estrogen receptor alpha (ESR1) in Chinese alligator, Alligator sinensis. Comp Biochem Physiol B Biochem Mol Biol. 2016;200:28–35.
Zhang R, Zhang SZ, Zhu X, Zhou YK, Wu XB. Molecular characterization of the Chinese alligator follicle-stimulating hormone beta subunit (FSH beta) and its expression during the female reproductive cycle. Comp Biochem Phys B. 2015;183:49–57.
Zhang R, Zhang S, Zhu X, Zhou Y, Wu X. Follicle-stimulating hormone receptor (FSHR) in Chinese alligator, Alligator sinensis: molecular characterization, tissue distribution and mRNA expression changes during the female reproductive cycle. Anim Reprod Sci. 2015;156:40–50.
Zhang R, Nie H, Duan S, Yan P, Izaz A, Wang R, Zhou Y, Wu X. Cloning, characterization and expression profile of kisspeptin1 and the kisspeptin1 receptor at the hypothalamic-pituitary-ovarian axis (HPO) of Chinese alligator, Alligator sinensis, during the reproductive cycle. Reprod Fertil Dev. 2020;32:792–804.
Sun HJ, Zuo XB, Sun L, Yan P, Zhang F, Xue H, Li E, Zhou YK, Wu R, Wu XB. Insights into the seasonal adaptive mechanisms of Chinese alligators (Alligator sinensis) from transcriptomic analyses. Aust J Zool. 2018;66(2):93–102.
Cooper ST, Sell SS, Fahrenkrog M, Wilkinson K, Howard DR, Bergen H, Cruz E, Cash SE, Andrews MT, Hampton M. Effects of hibernation on bone marrow transcriptome in thirteen-lined ground squirrels. Physiol Genomics. 2016;48(7):513–25.
Faherty SL, Villanueva-Canas JL, Klopfer PH, Alba MM, Yoder AD. Gene expression profiling in the hibernating primate, Cheirogaleus medius. Genome Biol Evol. 2016;8(8):2413–26.
Hampton M, Melvin RG, Andrews MT. Transcriptomic analysis of brown adipose tissue across the physiological extremes of natural hibernation. PLoS One. 2013;8(12):e85157.
Lei M, Dong D, Mu S, Pan YH, Zhang S. Comparison of brain transcriptome of the greater horseshoe bats (Rhinolophus ferrumequinum) in active and torpid episodes. PLoS One. 2014;9(9):e107746.
Luan Y, Ou J, Kunze VP, Qiao F, Wang Y, Wei L, Li W, Xie Z. Integrated transcriptomic and metabolomic analysis reveals adaptive changes of hibernating retinas. J Cell Physiol. 2018;233(2):1434–45.
Nespolo RF, Gaitan-Espitia JD, Quintero-Galvis JF, Fernandez FV, Silva AX, Molina C, Storey KB, Bozinovic F. A functional transcriptomic analysis in the relict marsupial Dromiciops gliroides reveals adaptive regulation of protective functions during hibernation. Mol Ecol. 2018;27(22):4489–500.
Wan QH, Pan SK, Hu L, Zhu Y, Xu PW, Xia JQ, Chen H, He GY, He J, Ni XW, et al. Genome analysis and signature discovery for diving and sensory properties of the endangered Chinese alligator. Cell Res. 2013;23(9):1091–105.
Davis GH. Fecundity genes in sheep. Anim Reprod Sci. 2004;82-83:247–53.
Rothschild MF, Messer L, Day A, Wales R, Short T, Southwood O, Plastow G. Investigation of the retinol-binding protein 4 (RBP4) gene as a candidate gene for increased litter size in pigs. Mamm Genome. 2000;11(1):75–7.
Draper BW, McCallum CM. Moens CB: nanos1 is required to maintain oocyte production in adult zebrafish. Dev Biol. 2007;305(2):589–98.
Shindo T, Kurihara H, Kuno K, Yokoyama H, Wada T, Kurihara Y, Imai T, Wang Y, Ogata M, Nishimatsu H, et al. ADAMTS-1: a metalloproteinase-disintegrin essential for normal growth, fertility, and organ morphology and function. J Clin Invest. 2000;105(10):1345–52.
Bole-Feysot C, Goffin V, Edery M, Binart N, Kelly PA. Prolactin (PRL) and its receptor: actions, signal transduction pathways and phenotypes observed in PRL receptor knockout mice. Endocr Rev. 1998;19(3):225–68.
Kashimada K, Pelosi E, Chen H, Schlessinger D, Wilhelm D, Koopman P. FOXL2 and BMP2 act cooperatively to regulate follistatin gene expression during ovarian development. Endocrinology. 2011;152(1):272–80.
Lin JQ, Zhou Q, Yang HQ, Fang LM, Tang KY, Sun L, Wan QH, Fang SG. Molecular mechanism of temperature-dependent sex determination and differentiation in Chinese alligator revealed by developmental transcriptome profiling. Sci Bull. 2018;63(4):209–12.
Shoemaker CM, Crews D. Analyzing the coordinated gene network underlying temperature-dependent sex determination in reptiles. Semin Cell Dev Biol. 2009;20(3):293–303.
Jorgensen CB. Relations between hibernation and ovarian functions in a temperate zone frog, Rana-Temporaria. Acta Zool. 1984;65(4):239–47.
Kuang Z, Yao Y, Shi Y, Gu Z, Sun Z, Tso J. Winter hibernation and UCHL1-p34cdc2 association in toad oocyte maturation competence. PLoS One. 2013;8(10):e78785.
Hossain MM, Sohel MM, Schellander K, Tesfaye D. Characterization and importance of microRNAs in mammalian gonadal functions. Cell Tissue Res. 2012;349(3):679–90.
Luense LJ, Carletti MZ, Christenson LK. Role of dicer in female fertility. Trends Endocrinol Metab. 2009;20(6):265–72.
Liu H-C, Tang Y, He Z, Rosenwaks Z. Dicer is a key player in oocyte maturation. J Assist Reprod Genet. 2010;27(9–10):571–80.
Suh N, Baehner L, Moltzahn F, Melton C, Shenoy A, Chen J, Blelloch R. MicroRNA function is globally suppressed in mouse oocytes and early embryos. Curr Biol. 2010;20(3):271–7.
Sirotkin AV, Ovcharenko D, Grossmann R, Laukova M, Mlyncek M. Identification of microRNAs controlling human ovarian cell steroidogenesis via a genome-scale screen. J Cell Physiol. 2009;219(2):415–20.
Kucharski R, Maleszka J, Foret S, Maleszka R. Nutritional control of reproductive status in honeybees via DNA methylation. Science. 2008;319(5871):1827–30.
Amarasinghe HE, Clayton CI, Mallon EB. Methylation and worker reproduction in the bumble-bee (Bombus terrestris). Proc R Soc B-Biol Sci. 2014;281(1780):20132502.
Si Y, Ding Y, He F, Wen H, Li J, Zhao J, Huang Z. DNA methylation level of cyp19a1a and Foxl2 gene related to their expression patterns and reproduction traits during ovary development stages of Japanese flounder (Paralichthys olivaceus). Gene. 2016;575(2):321–30.
Ding Y, He F, Wen H, Li J, Ni M, Chi M, Qian K, Bu Y, Zhang D, Si Y. DNA methylation status of cyp17-II gene correlated with its expression pattern and reproductive endocrinology during ovarian development stages of Japanese flounder (Paralichthys olivaceus). Gene. 2013;527(1):82–8.
Bai J, Gong W, Wang C, Gao Y, Hong W, Chen SX. Dynamic methylation pattern of cyp19a1a core promoter during zebrafish ovarian folliculogenesis. Fish Physiol Biochem. 2016;42(3):947–54.
Zhang Y, Li F, Feng X, Yang H, Zhu A, Pang J, Han L, Zhang T, Yao X, Wang F. Genome-wide analysis of DNA methylation profiles on sheep ovaries associated with prolificacy using whole-genome bisulfite sequencing. BMC Genomics. 2017;18(1):759.
Calicchio R, Doridot L, Miralles F, Mehats C, Vaiman D. DNA methylation, an epigenetic mode of gene expression regulation in reproductive science. Curr Pharm Des. 2014;20(11):1726–50.
Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.
Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.
Anders S, Pyl PT, Huber W. HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166–9.
Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26:136–8.
Chen C, Chen H, Zhang Y, Thomas HR, Frank MH, He Y, Xia R. TBtools: an integrative toolkit developed for interactive analyses of big biological data. Mol Plant. 2020;13(8):1194–202.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.
Kalvari I, Argasinska J, Quinones-Olvera N, Nawrocki EP, Rivas E, Eddy SR, Bateman A, Finn RD, Petrov AI. Rfam 13.0: shifting to a genome-centric resource for non-coding RNA families. Nucleic Acids Res. 2018;46(D1):D335–42.
Smit AFA, Hubley R, Green P. RepeatMasker Open-4.0. 2013-2015. http://www.repeatmasker.org.
Wen M, Shen Y, Shi S. Tang T: miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformatics. 2012;13:140.
Friedlander MR, Mackowiak SD, Li N, Chen W. Rajewsky N: miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40(1):37–52.
Betel D, Koppal A, Agius P, Sander C, Leslie C. Comprehensive modeling of microRNA targets predicts functional non-conserved and non-canonical sites. Genome Biol. 2010;11(8):R90.
Kertesz M, Iovino N, Unnerstall U, Gaul U, Segal E. The role of site accessibility in microRNA target recognition. Nat Genet. 2007;39(10):1278–84.
Rehmsmeier M, Steffen P, Hochsmann M, Giegerich R. Fast and effective prediction of microRNA/target duplexes. RNA. 2004;10(10):1507–17.
Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for bisulfite-Seq applications. Bioinformatics. 2011;27(11):1571–2.
Wang Z, Li X, Jiang Y, Shao Q, Liu Q, Chen B. Huang D: swDMR: a sliding window approach to identify differentially methylated regions based on whole genome bisulfite sequencing. PLoS One. 2015;10(7):e0132866.
Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14.
Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21(19):3787–93.
We would like to thank Li-Ming Fang, Zhen-Wei Wang, Wei-Qiang Zou, Da-Bin Ren and Ju-Min Xu (Changxing Yinjiabian Chinese Alligator Nature Reserve) for their help in sample collection.
This work was supported by the Key Program of the National Natural Science Foundation of China (31530087), the National Key Program (2016YFC0503200) from the Ministry of Science and Technology of China, a special grant from the State Forestry Administration, and the Fundamental Research Funds for the Central Universities of China.
Ethics approval and consent to participate
All the experimental procedures with Chinese alligator in this study had been given prior approval by the Animal Ethics Committee of Zhejiang University (ZJU2015–154-13) and the State Forestry Administration of China [Forest Conservation Permission Document (2014) 1545].
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Most enriched KEGG pathways in DEGs between SF_OVA_R and other ovary samples.
S1 Expression patterns of genes involved in sex differentiation and fertility in ovary samples from Chinese alligators overwintering in different environments. a. Expression alterations of genes involved in sex differentiation and fertility in the gonads of female alligators overwintering in the warm room. b. Expression heatmap of DEGs involved in sex differentiation and fertility in the ovary samples. c. Expression patterns of GDF9, PRLRB, BMP2B, SOX9, StAR, and FGFR2A in the ovary samples.
Correlations between DEGs and DMGs in ovaries in the breeding season. a. Venn diagram of DEGs and DMGs. b. Venn diagram of up- and downregulated DEGs and hyper- and hypomethylated DMGs. c, d. Venn diagrams of up- and downregulated DEGs and hypomethylated DMGs with DMR in promoters (c) and gene bodies (d). e, f. Venn diagrams of up- and downregulated DEGs and hypermethylated DMGs with DMR in promoters (e) and gene bodies (f).
DMGs and/or DEGs between SF_OVA_R and SF_OVA samples in sex determination, fertility, and oocyte meiosis pathways.
The most enriched (FDR < 0.1) pathways of season-biased DEGs in the testis of Chinese alligator.
The most enriched (FDR < 0.05) GO terms of season-biased DEGs in the testis of Chinese alligator.
Gene expression related to the androgen biosynthesis pathway in the testes. a. Androgen biosynthesis pathway. Protein names in orange ovals indicate enzymes encoded by genes that are upregulated in testes collected in summer. b. Expression levels of AMH, MTNR1A, CYP11A, HSD3B, and CYP17 in the testes.
Regulation of DEGs and DMGs encoding transcription factors (TFs). a. Venn diagram of TFs and DMGs. b. Venn diagram of TFs, DEGs, and DMGs. c. Associations between TFs and key genes with ovary-related functions.
About this article
Cite this article
Lin, JQ., Yu, J., Jiang, H. et al. Multi-omics analysis reveals that natural hibernation is crucial for oocyte maturation in the female Chinese alligator. BMC Genomics 21, 774 (2020). https://doi.org/10.1186/s12864-020-07187-5
- DNA methylation
- Hibernation environment