Heterosis of growth trait regulated by DNA methylation and miRNA in allotriploid fish
Epigenetics & Chromatin volume 15, Article number: 19 (2022)
Heterosis of growth traits in allotriploid fish has benefited the production of aquaculture for many years, yet its genetic and molecular basis has remained obscure. Now, an allotriploid complex, including two triploids and their diploid inbred parents, has provided an excellent model for investigating the potential regulatory mechanisms of heterosis.
Here, we performed a series of analyses on DNA methylation modification and miRNA expression in combination with gene expression in the allotriploid complex. We first established a model of cis- and trans-regulation related to DNA methylation and miRNA in allotriploids. Then, comparative analyses showed that DNA methylation contributed to the emergence of a dosage compensation effect, which reduced gene expression levels in the triploid to the diploid state. We detected 31 genes regulated by DNA methylation in the subgenomes of the allotriploids. Finally, the patterns of coevolution between small RNAs and their homoeologous targets were classified and used to predict the regulation of miRNA expression in the allotriploids.
Our results uncovered the regulatory network between DNA methylation and miRNAs in allotriploids, which not only helps us understand the regulatory mechanisms of heterosis of growth traits but also benefits the study and application of epigenetics in aquaculture.
Allotriploids, including two or more different genomes in somatic cells, always exhibit phenotypic differences in growth [1, 2], gonadogenesis , and innate immunity [4, 5] from diploid relatives or parents. Some phenotypic diversities, such as heterosis, provide us with a useful research resource for investigating their potential mechanisms and could improve agricultural production by engineering breeding [6, 7]. Some potential molecular regulatory mechanisms involving genomic  and transcriptomic changes [9, 10] have been reported as key factors in shaping heterosis. However, studies on epigenetic modification and post-transcriptional modification are very rare in allotriploids. How modifications of homoeologs affect phenotypes in allotriploids is interesting but largely unknown.
Carassius auratus red var. and Cyprinus carpio L., belonging to different genera of the subfamily Cyprinidae, experienced a common whole-genome duplication event approximately 13.75 million years ago (Mya) . Their divergence time was dated to 10.0 Mya . As important aquaculture species in China, the annual production capacity of common carp and crucian carp in 2018 was approximately 2.9 and 2.7 million tons, respectively . Common carp and crucian carp are mainly distributed in Asia and Europe, reaching as far as the Arctic Circle. Their distinct morphological characteristics, including body color and shape [13, 14], contribute to their extensive adaptability to a diversified environment of slow-moving rivers, lakes, and ponds. Heterosis of the growth trait was observed in the allotriploid complex of common carp and crucian carp [7, 9]. Among the two triploids, 3nRC2 had two sets of subgenome C (originating from common carp) and one set of subgenome R (originating from red crucian carp) [7, 9], while 3nR2C had one set of subgenome C and two sets of subgenome R [7, 9].
Dosage compensation is always detected in the two X chromosomes of female humans, leading to the silencing of one of the X chromosomes . Moreover, some studies have reported dosage compensation in allotriploid animals and plants, reducing transcript levels to the diploid state [1, 16, 17]. Incomplete dosage compensation and maternal effects were predicted in the two triploid fish (3nR2C and 3nRC2) based on the gene expression changes between the inbred parents and the two triploids . This also had various effects on the expression of growth-regulated genes, resulting in heterosis of the growth trait. However, the effects of epigenetic regulation on these expression changes are still unknown. Dosage compensation on the sex chromosomes is not only regulated by DNA methylation  but is also affected by microRNAs (miRNAs), noncoding RNAs that regulate gene expression at the post-transcriptional level [19, 20]. However, the expansion of miRNA families and their target sequences, arising through gene duplication, complicates the miRNA–mRNA regulatory network in allotriploids. Although some studies have reported miRNA expression changes between triploid and diploid individuals, the expression of miRNAs and the target genes did not distinguish between different homoeologs in subgenomes . The two homoeologs in 3nR2C and 3nRC2 originate from species of different genera and exhibit nucleotide differences in parts of the coding sequences and their regulatory regions. Therefore, the specific recognition and detection of mRNA and miRNA expression levels and DNA methylation loci will improve our understanding of the regulatory interactions between the two subgenomes in intergeneric hybrids.
The specific genetic status of the two allotriploids (3nR2C and 3nRC2) offers us a new insight into how dosage compensation occurs between two sets of the same chromosomes [9, 16, 17], and helps us further understand the relationship between genetic changes and heterosis of growth traits in allotriploids. Here, we performed DNA methylation sequencing (methyl-seq), miRNA-seq, and mRNA-seq in the muscles of the two triploids and their inbred parents (C. auratus red var. and C. carpio L.). Integrating analyses of gene expression with DNA methylation and miRNA expression, the changes between diploids and triploids or between subgenomes of 3nR2C and 3nRC2 will help us investigate potential mechanisms of growth traits in the triploids and improve the aquaculture of common carp and crucian carp by epigenetic engineering and breeding.
Heterosis of growth traits in allotriploid fish
In our studies, the two triploids (3nR2C and 3nRC2) were obtained by backcrossing female allotetraploids of C. auratus red var. and C. carpio L. (4nR2C2) with male C. auratus red var. (2nRR) and C. carpio L. (2nCC), respectively (Fig. 1A) [7, 9]. Eight months after hatching, some growth traits, including body length (BL), body height (BH), height of back muscle (HBM), and body weight (BW), showed higher growth rates in the two triploids (3nR2C and 3nRC2) than in the diploid inbred parents (2nRR and 2nCC) (Fig. 1). Among these traits, similar values in HBM were observed between the two triploids (Fig. 1). Moreover, the ratios of BL vs. BH in the two triploids exhibited intermediate characteristics between the inbred parents (Fig. 1). Specifically, the two triploids exhibited a predilection for their maternal parents for their distinct growth phenotypes involving HBM and ratios of BL vs. BH (Fig. 1A, B). These differences may be related to their genotypes and play a key role in heterosis.
Cis- and trans-regulation models of DNA methylation modification and miRNA regulation
To investigate potential genetic mechanisms and their roles in the heterosis of growth trait heterosis, methyl-seq, miRNA-seq, and mRNA-seq data were obtained from 3nR2C, 3nRC2, and the two inbred parents. The complex regulation of gene expression in the two allotriploid fish led us to establish a model for investigating cis- and trans-regulation in DNA methylation modification and trans-regulation in miRNA regulation. In the allotriploid system, the expression of homoeologs R (originating from 2nRR) and C (originating from 2nCC) was not just regulated by their respective methylation modifications but also by similar miRNA systems in the two subgenomes (Fig. 2). The development of diagrammatic regulatory systems aided us in better understanding cis- and trans-regulation in allotriploids, as well as their effects on phenotype. Cis- and trans-regulation in hybrids was assessed as described by Wang et al. [22, 23]. The changes in DNA methylation levels and the expression of their regulated genes from the inbred parents to their triploids are used to assess DNA methylation-mediated cis-regulation. MicroRNA-mediated trans-regulation in hybrids is then regulated by trans-acting factors involving microRNA, which was assessed by comparing the change trends in microRNA expression and targeted gene expression between inbred parents and their triploids. The combined genomes of 2nRR and 2nCC were used as a reference genome for the two triploid fish in our subsequent analyses for three reasons: (1) the lack of an assembly genome for the triploids; (2) completed sets of original chromosomes from the two inbred parents ; and (3) research on conserved homoeologous genes in subgenomes. Therefore, analyses of orthologous gene pairs between 2nRR and 2nCC, and homoeologous gene pairs between subgenomes R and C in the two triploids (3nR2C and 3nRC2), were performed to study the interactions between DNA methylation and miRNA and mRNA expression in the triploids.
Methylation divergence between subgenomes R and C
To uncover DNA methylation differences from diploids to allotriploids, we performed whole-genome bisulfite sequencing (WGBS) in the muscles of 2nRR, 2nCC, 3nR2C, and 3nRC2 fish. After mapping the reads of the DNA methylation data to the 2nRR and 2nCC combined genome, 270.85 Gb (48.91%) of data uniquely mapped in the two triploids (Additional file 1: Table S1). Based on the mapped results, 138.39 Gb uniquely mapped data (49.33% of all methylation data) were detected in the three 2nCC samples, while 150.85 Gb data (54.73%) were found in 2nRR (Additional file 1: Table S1). Focusing on cytosines, the methylation levels in CpG islands (70.72%) were higher than those in CHH (0.30%) and CHG (0.30%) across all samples. However, the levels of methylation divergence at CpG and nonCpG sites were conserved among the two triploids and the inbred parents (Additional file 1: Table S2).
Although similar genome sizes were detected between the 2nRR (1.42 Gb) and 2nCC (1.49 Gb) genomes, different distributions of transposable elements (TEs) between them reflected that TEs had contributed to the diversified species divergence after a common allotetraploidization in their ancestors [11, 24] (Additional file 1: Table S3 and Fig. S1A). More TEs were found in 2nRR (38.99%) than in 2nCC (34.76%). The top three TE difference ratios between 2nRR and 2nCC were observed in orthologous chromosome pair (OCP) 43 (16.75%), 2 (10.04%), and 5 (10.06%). Meanwhile, high DNA methylation was observed in the region of the TEs (Additional file 1: Fig. S1B and C). We first identified global TE methylation levels and found that the TE methylation levels in 3nR2C trended toward 2nRR, while those in 3nRC2 trended toward 2nCC (Additional file 1: Fig. S1B). Then, we further detected the TE methylation levels in different subgenomes in the triploids and the respective ones in the inbred parents. Higher TE DNA methylation was observed in both subgenomes of 3nR2C than in 3nRC2 (Additional file 1: Fig. S1B).
Genome-wide analyses showed that methylation levels were higher in 2nRR than in 2nCC (p < 0.01 in t-test, Fig. 3A and Additional file 1: Table S4). These results showed that DNA methylation divergence occurred between 2nRR and 2nCC after their shared ancestral allotetraploidization [11, 24]. Meanwhile, these DNA methylation divergences in the inbred parents were also observed between subgenomes R and C in the two triploids (Fig. 3A and Additional file 1: Table S4). These results suggest high heritability for DNA methylation. Then, the DNA methylation differences between the four fish were mainly enriched in the first exon coding region, 2 kb upstream of the transcription start site (TSS), and 2 kb downstream of the transcription termination site (TTS) (Additional file 1: Fig. S2). Moreover, clustering analysis in the subgenome showed that subgenome R in the two triploids clustered with 2nRR, yet subgenome C clustered with 2nCC (Additional file 1: Fig. S2). Similar to the TE methylation levels, the methylation levels of 2 kb upstream of TSS (Up 2 k) and 2 kb downstream of TTS (Down 2 k) were higher in the two triploids than in the inbred parents (p < 0.01), except those in subgenome C of 3nRC2, in which the DNA methylation levels were close to those of the parental 2nCC (Fig. 3A). These results shed light on high DNA methylation and low transcriptional efficiency in triploids.
Incomplete dosage compensation led by DNA methylation
Recent studies have reported incomplete dosage compensation . To investigate whether DNA methylation is a potential cause of this, we performed integrated analyses of differential methylation and gene expression. First, genotype analysis showed two sets of 2nRR chromosomes and one set of 2nCC chromosomes in 3nR2C (genotype type: AAB), and two sets of 2nCC chromosomes and one set of 2nRR chromosomes in 3nRC2 (genotype type: ABB) . The two sets of R genomes/subgenomes were in 2nRR and 3nR2C, whereas the two sets of C genomes/subgenomes were in 2nCC and 3nRC2. However, the gene expression of R homoeologs in 3nR2C (× 1 in 3nR2C) was significantly lower than that in 2nRR (× 1) (p < 0.01 in t-test, Fig. 3B). Based on the same regulatory systems, we doubled the gene expression values of the R and C subgenomes in the two allotriploids (3nR2C and 3nRC2) (× 2). Then, we found that the doubled values of the R subgenome in 3nR2C (× 2 in 3nR2C) were equal to those in 2nRR (× 1) (p > 0.05, Fig. 3B). A similar phenomenon was observed between the doubled values of the C subgenome in 3nR2C (× 2 in 3nRC2) and those in 2nCC (× 1) (p > 0.05, Fig. 3B). So, we inferred that the four sets of R and C subgenomes in 3nR2C and 3nRC2 can transcribe the same dosage of mRNAs as the parental 2nRR and 2nCC based on the regulatory models in the two triploids.
In analyses of differentially methylated regions (DMRs) between the two triploids and their inbred parents, most DMRs in subgenome R were detected in the comparison “3nRC2 vs. 2nRR” (11,152), while the fewest DMRs were observed in the comparison “3nR2C vs. 2nRR” (Fig. 3C). In subgenome C, most DMRs were observed in the comparison “3nR2C vs. 3nRC2” (3,999), while the fewest DMRs were observed in the comparison “3nRC2 vs. CC” (Fig. 3C). Obviously, many DMRs were distributed between the unequal subgenomes of R and C. In the comparison “3nR2C vs. 3nRC2,” we observed that more gene copies always led to higher methylation levels in the triploids (Fig. 3C, D). Specifically, these effects were more obvious in subgenome R than in subgenome C (Fig. 3C). In the comparisons of triploids and diploids, higher methylation levels and a greater number of DMGs reflected that whole genome-wide methylation modification inhibited transcription of both homoeologs R and C and led to the triploid expression levels being close to the diploid states.
DNA methylation changes accompanied by homoeologous recombination
We identified 165 hypermethylated (DMGsA–B > 0.6) and 740 hypomethylated genes (0.4 < DMGsA–B < 0.6) based on the methylation differences between 2nRR and 2nCC. Among these genes, 15 hyper-DMGs and 114 hypo-DMGs were shared in both the triploids, suggesting that the potential DNA methylation changes may have originated from the parental 4nR2C2 (Fig. 3E). These DMGs were considered potential homoeologous recombinant genes (HRGs) because homoeologous recombination (HR) always leads to a decrease in DNA methylation divergence between homoeologs R and C in promoter regions (Fig. 3E). For assessing the accuracy of HRGs predicted by DNA methylation, the mapped read number of homoeologs R and C was obtained from whole-genome resequencing, and unequal HRs were detected in some genes (52 in 3nR2C, 50 in 3nRC2, red dots, Fig. 3F) based on the ratios of the homoeologs R vs. C . Meanwhile, a similar shape (blue line, Fig. 3F) between 3nR2C and 3nRC2 indicated a similar distribution of HRs in these genes, suggesting that the prediction that the HRs were inherited from paternal 4nR2C2 was highly accurate.
Gene expression regulated by DNA methylation
DNA methylation in the promoter region is a major epigenetic modification and is always negatively correlated with gene expression. In the hybrid system of the two triploids, 17 and 14 genes showed negative correlation between differential expression (DE) and differential methylation (DM) in subgenomes R and C, respectively (Fig. 4A–B). These genes contributed to multiple functions, including regulation of response to stimulus and cellular metabolic process (Additional file 1: Fig. S3). Moreover, comparing the two triploids with their inbred diploid parents, various genes were identified with negative correlations between the DE and DM values (red dots, Additional file 1: Fig. S4), while other genes showed a positive correlation (black dots, Additional file 1: Fig. S4). Integrated analyses of methylation and gene expression changes demonstrated that only a few genes were regulated by DNA methylation, while most of the expressed genes might be regulated by trans-regulation in the hybrid systems.
Patterns of miRNA regulation
In a hybrid system, miRNA regulation is a key trans-regulatory mechanism that may play a dominant role in gene expression regulation. We obtained 129.09 million clean reads (2.91 Gb) from small RNA-seq from the four fish with three biological replicates (Additional file 1: Table S5). Among these data, 121.10 million reads mapped to the combined genomes of the inbred parents, and 258 miRNAs were detected in 2nRR, 2nCC, and the two triploids using the zebrafish miRNA. The multiple transcription sites of some miRNAs were not just detected in 2nRR and/or 2nCC but also in the subgenomes of R and C in the two triploids. Under these conditions, we established six patterns for classifying these conserved miRNAs, which might result from random duplication of chromosome regions or whole-genome duplication (Fig. 5A). For the allotriploids, miRNAs in Pattern 1 (114) and Pattern 4 (51) possessed four and two copies, respectively, that were evenly distributed in subgenomes R and C (Fig. 5A). In addition, three copies were observed in 65 (Pattern 2) and 34 (Pattern 3) miRNAs, in which two miRNA copies were distributed in subgenome C or R, respectively, and the other copy was distributed in the other subgenome (Fig. 5A). The miRNAs in Pattern 5 had two copies distributed in only subgenome R (26) or subgenome C (6) (Fig. 5A). Pattern 6 describes miRNAs with only a single copy in subgenome R (133) or C (69) (Fig. 5A). The evolutionarily conserved class of miRNAs sheds light on the relationship between miRNA copy number and their expression.
Because of the evolutionarily conserved miRNAs, the expression of each miRNA was calculated by combining the expression values of subgenomes R and C based on the mapped files and annotation results. The analysis of differential expression (DE) was performed on 256 miRNAs, of which 33 miRNAs were detected between 3nRC2 and 2nCC, and 62 miRNAs were detected between 3nR2C and 2nCC (Fig. 5B, C and Additional file 1: Fig. S5A). DE was detected in 28 miRNAs between 3nRC2 and 2nRR, and for only two miRNAs between 3nR2C and 2nRR (Fig. 5D and Additional file 1: Fig. S5A). Meanwhile, 26 miRNAs showed DE between 3nRC2 and 3nR2C (Fig. 5E). Among these miRNAs, dre-mir-182 and dre-mir-183 showed DE in the three comparisons (Additional file 1: Fig. S5B).
Contribution of DNA methylation modification and miRNA regulation to heterosis of growth traits
Our aims were to investigate the genetic mechanisms of DNA methylation and miRNA in triploids and to further understand their potential effects on heterosis of growth traits. The functional distribution of the hyper-/hypo-DMGs was performed and annotated with Gene Ontology databases and 97 pathways (Additional file 1: Table S6 and Fig. S6). Among these genes, six genes were enriched for the growth hormone synthesis, secretion, and action pathway (map04935), while nine genes were distributed in the insulin secretion pathway (map04911) (Additional file 1: Figs. S7, 8). Meanwhile, high methylation levels led to the expression silencing of homoeolog C in ATF4 and SNAP25 of the two triploids, while low methylation levels and expression of homoeolog R were detected (Fig. 4E). On the contrary, low methylation levels of homoeolog R were observed in CREB1 and STAT5A, while high methylation levels were detected in homoeolog C (Fig. 4E). Some growth-regulated genes were regulated by miRNA regulation, especially in the growth hormone synthesis, secretion, and action pathway (map04935). Among these growth-regulated genes, the expression levels of P85B and EP300 were regulated by mir-205, which showed DE between the two triploids (Fig. 5E). However, no DE of the two genes between the two triploids indicated a complex regulation system in the triploids. The heterosis of growth traits in the triploids might involve cooperative changes to a series of growth-regulated genes.
Both common carp and crucian carp are important aquaculture species in China because of their wide distribution and adaptation [12,13,14]. The two triploid fish, which have been obtained via the interploid crossing of 4nR2C2 with its inbred parents (common carp and red crucian carp), are also important aquaculture varieties in the Yangtze River because of their stronger innate immunity and faster growth rate than their inbred parents (Fig. 1) [7, 9]. To investigate the potential regulatory mechanisms in the triploids and their effects on the growth traits, we focused on cis- and trans-regulation (DNA methylation modification and miRNA regulation) in combination with gene expression (Fig. 2).
The mRNAs and miRNAs had more gene copies in the triploids than in the diploid inbred parents, while heterosis of the growth rate and body size was observed in the triploids (Fig. 1). However, a dosage compensation effect caused the gene expression of the three sets of subgenomes in the triploids to be close to those of the two sets of subgenomes in the inbred parents (Fig. 3). This was regulated by DNA methylation and other regulatory mechanisms, and effectively balanced the steady state of cells and individuals [9, 17]. However, the gene expression changes caused by the effect were not identical. Higher expression of GH and IGF1 in the triploid progenies compared with that in the diploid parents may contribute to the heterosis of growth . Meanwhile, the coordinated regulation of DNA methylation and miRNA expression contributed to various effects on the expression of targeted growth-regulated genes (Figs. 2–5), which may be related to the heterosis of growth in the triploid fish.
Specific genotypes (two sets of subgenomes R and one set of subgenome C in 3nR2C, one set of subgenome R and two sets of subgenomes C in 3nRC2) and different growth phenotypes also shed light on the complex regulatory mechanisms of heterosis of growth traits. Functional divergences between the two inbred parents (such as growth rate and body size, 2nCC > 2nRR) revealed their divergent regulatory systems (Figs. 1–5). In the hybrid system of the triploids, the expression of targeted mRNAs was regulated by miRNAs originating from both subgenomes R and C [16, 27, 28]. Comparative analysis of the two allotriploids revealed that more gene copies of subgenome C may contribute to the expression bias of homoeolog C and further benefit the large body size and rapid growth rate (Fig. 1). Under these conditions, exchanges of the promoter region between subgenomes R and C could disturb the DNA methylation originating from the two inbred parents and further shape the diversified expression of homoeologs (Fig. 3) [25, 29]. However, the dispersed distribution of these recombination events in chromosome segments may result from DNA damage response and repair (via recombination) . Additional studies on the recombination mechanisms in the allotetraploid fish are necessary. Meanwhile, various changes in miRNA expression from the inbred parents to the allotriploid progeny also contributed to the complexity of homoeolog expression (Fig. 5).
In the triploids, gene expression changes from the inbred parents and expression bias of homoeologs play a key role in heterosis of growth . Meanwhile, the different heterosis of growth traits may be relative to the different dosage of R and C subgenomes in the two triploids (2:1 in 3nR2C and 1:2 in 3nRC2), which derive from the distinctly different body sizes and ratios of BL vs. BH between 2nRR and 2nCC. A significant effect of different subgenomes on heterosis was also described in triploid maize . These expression changes were regulated by complex cis- and trans-regulatory mechanisms, including DNA methylation and miRNAs. Our results help us to further understand the potential regulatory mechanisms of heterosis of growth traits in triploid fish.
The fish used in this study were diploid C. auratus red var. (2nRR), diploid C. carpio L. (2nCC), an triploid fish (3nR2C) obtained from hybridization between C. auratus red var. (♀) and an allotetraploid of C. auratus red var. x C. carpio L. (4nR2C2, ♂), and an triploid fish (3nRC2) obtained from hybridization between C. carpio (♀) and 4nR2C2 (♂). The fish were fed in a pool located in the Engineering Center of Polyploidy Fish Breeding of the National Education Ministry, Hunan Normal University, China. Two-year-old male individuals of each fish were collected with three biological replicates. These fish were deeply anesthetized with 300 mg/L Tricaine Methanesulfonate (Sigma-Aldrich, St Louis, MO, USA) for 10 min (25 °C) in a separate tank. After confirming death, they were collected from the water for dissection. A flow cytometer was used to measure the DNA content of the erythrocytes. To avoid contamination, all tissue samples were collected using sterilized scissors and tweezers.
Identification of orthologous chromosomes and orthologous genes
Despite the lack of a reference genome for the nascent triploids, the combined genome of the two inbred parents could still be used as the in silico genome because of the presence of the two subgenomes of the inbred parents in the triploids. Therefore, the combined genome was used as the reference genome for investigating DNA methylation and microRNA (miRNA) and mRNA expression. Then, orthologs between C. auratus red var. and C. carpio L. were considered the homoeologs of subgenomes in the two triploids (3nR2C and 3nRC2). Orthologous chromosomes between C. auratus red var. and C. carpio L. were detected using MCScanX with the thresholds of e-value < 10−5 and max gaps < 25 . We performed reciprocal highest scores in BLAST similarity searches of all annotated genes of C. auratus red var. and C. carpio L. with the criteria of e-value < 10−5, a minimum of 80% sequence length coverage, and 300 bp length of sequence identity. Then, high confidence orthologous gene pairs were finally identified based on the shared gene pairs satisfying the results of the reciprocal BLAST and orthologous chromosomes.
Whole-genome bisulfite sequencing and mapping of methyl-seq reads
The WGBS libraries were constructed following the standard protocol . Briefly, genomic DNA spiked-in with 0.5% unmethylated lambda DNA (Fermentas) was first fragmented into 200–300 bp fragments using a Covaris S220 system. After being concentrated, the sheared DNA was end-repaired, dA-tailed, and ligated with pre-methylated TruSeq DNA adapters (Illumina). Bisulfite conversion was conducted with an EZ DNA Methylation™ Kit (Zymo Research). After the bisulfite conversion, the converted templates were PCR amplified and quantified with Qubit 2.0. Then, the length of the insert fragment in the libraries was detected by Agilent 2100, and the effective concentration (> 2 nM) was quantified by qPCR assay. The final quality-ensured libraries were sequenced on a NovaSeq 6000 Sequencing System with paired end reads (2 × 150 bp) according to the standard protocol.
After quality checking of the methyl-seq reads, the clean reads of C. auratus red var. and C. carpio L. were mapped to the respective genomes, and the clean reads of the two triploids were mapped to the combined genome sequences of C. auratus red var. and C. carpio L. (Yellow River carp) [11, 24]. The Bismark analysis pipeline was used to detect the methylated loci with the mapped parameters (–score_min L, 0, − 0.2 − X 1000 –no-mixed–no-discordant) [33, 34]. The clean reads were used for mapping to the reference genome four times, and only the reads that mapped to the same position of the reference genome each time were retained in our next analysis. A binomial distribution test was performed to identify 5-methylcytosine for each cytosine site. Then the potential methylation sites were checked with the thresholds of coverage > 4 × depth and false discovery rate (FDR) < 0.05.
Analysis of methylation level in transposons and differentially methylated regions (DMRs)
The average CpG methylation was detected in different gene regions, including upstream (a window size of 100 bp for 2 kb regions), the gene body, and downstream (a window size of 100 bp for 2 kb regions) of the coding regions. The transposons of 2nRR and 2nCC were predicted using RepeatMasker, and the average CpG methylation in the upstream and downstream transposon regions (2 kb) was calculated and plotted using R. The regions with different methylation were detected using MOABS . The R packages DSS and bsseq were used to call DMRs between homoeologous gene pairs in the two triploids and between orthologous genes in 2nRR and 2nCC based on the threshold of e-value < 10−5. The functional enrichment of DMRs was performed by annotating them with the Gene Ontology and Kyoto Encyclopedia of Genes and Genomes databases with a threshold of p-value < 0.05.
Association analyses of DNA methylation and homoeologous recombination
The promoter sequences (2 kb) of DMRs were identified for analyzing the differentially methylated genes (DMGs). The DMGs of orthologous gene pairs in the two inbred parents and homoeologous gene pairs in the hybrids were classified into the following two categories: (1) hyper-DMGs conforming to the thresholds of an absolute value of differences in the methylation ratio between 2nRR and 2nCC (|DMGsA–B|) > 0.6 and an absolute value of differences in the methylation ratios between the two homoeologs of the triploids (|DMGsAs–Bs|) < 0.3, and 2) hypo-DMGs conforming to the thresholds of 0 <|DMGsA–B|< 0.6 in the inbred parents and |DMGsAs–Bs|< 0.2 in the triploids.
miRNA identification and differential expression
After DNase treatment, total RNA (~ 10 µg) of each sample was quality checked using a 2100 Bioanalyzer based on the values of the A260/A280 and A260/A230 ratios. Small RNAs were obtained by 6% polyacrylamide gel electrophoresis on Novex™ TBE gels (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. After adding adaptors, RT-PCR with adaptor-specific primers was used to generate a small RNA library from the two inbred parents (2nRR and 2nCC) and their two triploids (3nR2C and 3nRC2). The small RNA library was then sequenced by Illumina HiSeq 2500 (Illumina, San Diego, CA, USA) according to the manufacturer’s protocol.
Adaptor removal and quality control were performed by Cutadapt (v. 1.7.1), fastx_toolkit (v. 0.0.13) and NGSQCToolkit (v. 2.3.2). A similar strategy to methyl-seq was used for the mapping of miRNA reads. Subsequently, clean reads of 15–26 bp in 2nRR, 2nCC, and the two triploids were mapped to the combined genomes of 2nRR and 2nCC using Bowtie (v. 1.1.2)  and miRDeep2 (https://anaconda.org/bioconda/mirdeep2, v. 184.108.40.206) . Then, the annotation of known miRNAs was performed by mapping to hairpin and mature miRNA sequences of zebrafish (http://www.mirbase.org/, Release 22.1). The annotated zebrafish miRNA files were downloaded from Ensembl database (http://dec2017.archive.ensembl.org/info/data/ftp/index.html). Then, the miRNA expression levels in the four fish were calculated using the script “quantifier.pl” in miRDeep2 software . The number of mapped reads was obtained from the output results. The analyses of DE were performed using edgeR package  based on thresholds of |log2 fold change|> 1 and FDR < 0.01 for three biological replicates. The target genes of the miRNAs were obtained in two ways: (1) collection from miRTarBase database (http://mirtarbase.cuhk.edu.cn/php/index.php), and (2) prediction from miranda software .
Copy number of miRNAs between red crucian carp and common carp
The miRNAs identified by the above analyses were used in our next analyses. Redundant sequences of 2nRR and 2nCC were deleted using CD-HIT (v. 4.8.1) software (https://github.com/weizhongli/cdhit/releases) with a threshold of sequence identity = 1.0. Then, the miRNA sequences were mapped to the combined genomes of 2nRR and 2nCC using the script “mapper.pl” of miRDeep2 software . Among the genome locations of the miRNAs, potential homologous miRNAs were considered those with a distribution between the orthologous chromosomes of 2nRR and 2nCC. Lastly, the multiple locations of these potential homologous miRNAs were used to investigate the copy number of the homologous miRNAs.
Expression of homoeologous genes
Muscle mRNA-seq data were downloaded from the Sequence Read Archive database in NCBI (accession numbers: SRS4475351–SRS4475354). These samples are the same as the ones in the miRNA and DNA methylation studies. Low-quality reads and adapters were trimmed off by Trimmomatic . Then, the quality of the Illumina reads was checked with FastQC software (v. 0.11.7) . A similar strategy to that of methyl-seq was used for mapping to investigate homoeologous gene expression. mRNA-seq reads were mapped to the combined genomes of C. auratus red var.  and C. carpio L.  using STAR software  with the parameters “–outSAMmapqUnique 255.” The unique mapped reads were collected using samtools  for the next expression analyses. Lastly, the read count was collected using HTSeq .
These genes were deleted in homoeolog expression analyses with the threshold of mapped read counts < 5 in each sample. The homoeolog expression of the two triploids was calculated based on the mapped read counts in the combined genome file, while expression levels of homologous genes in the two inbred parents were calculated based on the mapped results of the corresponding genome. The analyses of normalization and DE were performed using the edgeR package . Differentially expressed genes (DEGs) were detected with the thresholds of |log2 FC|> 1 and FDR < 5% with three biological replicates.
Obtaining mapped reads in homoeologous gene pairs using whole-genome resequencing
To investigate the homoeologous recombination (HRs) at the genomic level, whole-genome resequencing data were obtained from sequencing of HiSeq X Ten (paired end, 150 bp) in the NCBI Sequence Read Archive database . Quality checking and adapter removal were performed with Trim Galore (v. 0.4.0) (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) and Cutadapt (v. 1.2.1), respectively. Then, the high-quality reads were mapped to the combined genome of C. auratus and C. carpio using BWA (v. 0.7.17-r1188) with the default parameters. Then, the number of mapped reads in each gene was calculated using htseq-count (v. 0.12.4) with the threshold of “-m union–nonunique = none.” The HRs were assessed based on the ratio of uniquely mapped reads of homoeologs R vs. C .
Availability of data and materials
The raw reads of the DNA methylation data and small RNA sequencing data have been deposited in the National Genomics Data Center (NGDC) (https://bigd.big.ac.cn/gsub/) under the BioProject accession no. PRJCA003625. The mRNA sequencing and whole-genome resequencing data can be downloaded from NCBI Sequence Read Archive under the accession nos. SRR9983189, SRR9983190, SRR8735277–SRR8735279, SRR8712975–SRR8712976, SRR8712978–SRR8712991, SRR9203584, SRR9185089, and SRR9185090.
Million years ago
DNA methylation sequencing
Diploid C. auratus red var.
Diploid C. carpio L.
- 4nR2C2 :
An allotetraploid of C. auratus red var. x C. carpio L.
Triploid fish obtained from hybridization between C. auratus red var. (♀) and an allotetraploid of C. auratus red var. x C. carpio L. (4nR2C2, ♂)
- 3nRC2 :
Triploid fish obtained from hybridization between C. carpio (♀) and 4nR2C2(♂)
Whole-genome bisulfite sequencing
False discovery rate
Differentially methylated regions
Differentially methylated genes
Differentially expressed genes
Height of back muscle
Orthologous chromosome pair
Transcription start site
Transcription termination site
Homoeologous recombinant genes
National Genomics Data Center
Ren L, Tang C, Li W, Cui J, Tan X, et al. Determination of dosage compensation and comparison of gene expression in a triploid hybrid fish. BMC Genomics. 2017;18:38.
Song C, Liu S, Xiao J, He W, Zhou Y, et al. Polyploid organisms. Sci China. Life Sci. 2012;55:301–11.
Xu K, Wen M, Duan W, Ren L, Hu F, et al. Comparative analysis of testis transcriptomes from triploid and fertile diploid cyprinid fish. Biol Reprod. 2015;92:95.
Xiao J, Fu Y, Wu H, Chen X, Liu S, et al. MAVS of triploid hybrid of red crucian carp and allotetraploid possesses the improved antiviral activity compared with the counterparts of its parents. Fish Shellfish Immunol. 2019;89:18–26.
Ignatz EH, Braden LM, Benfey TJ, Caballero-Solares A, Hori TS, et al. Impact of rearing temperature on the innate antiviral immune response of growth hormone transgenic female triploid Atlantic salmon (Salmo salar). Fish Shellfish Immunol. 2020;97:656–68.
Huang X, Yang S, Gong J, Zhao Q, Feng Q, et al. Genomic architecture of heterosis for yield traits in rice. Nature. 2016;537:629–33.
Chen S, Wang J, Liu S, Qin Q, Xiao J, et al. Biological characteristics of an improved triploid crucian carp. Sci China C Life Sci. 2009;52:733–8.
Cassady JP, Young LD, Leymaster KA. Heterosis and recombination effects on pig reproductive traits. J Anim Sci. 2002;80:2303–15.
Ren L, Yan X, Cao L, Li J, Zhang X, et al. Combined effects of dosage compensation and incomplete dominance on gene expression in triploid cyprinids. DNA Res. 2019;26:485–94.
Rapp RA, Udall JA, Wendel JF. Genomic expression dominance in allopolyploids. BMC Biol. 2009;7:18.
Luo J, Chai J, Wen Y, Tao M, Lin G, et al. From asymmetrical to balanced genomic diversification during rediploidization: Subgenomic evolution in allotetraploid fish. Sci Adv. 2020;6:7677.
Shao GM, Li XY, Wang Y, Wang ZW, Li Z, et al. Whole genome incorporation and epigenetic stability in a newly synthetic allopolyploid of Gynogenetic Gibel Carp. Genome Biol Evol. 2018;10:2394–407.
Xu P, Zhang X, Wang X, Li J, Liu G, et al. Genome sequence and genetic diversity of the common carp, Cyprinus carpio. Nat Genet. 2014;46:1212–9.
Zhang Y, Liu J, Fu W, Xu W, Zhang H, et al. Comparative Transcriptome and DNA methylation analyses of the molecular mechanisms underlying skin color variations in Crucian carp (Carassius carassius L.). BMC Genet. 2017;18:95.
Nguyen DK, Disteche CM. Dosage compensation of the active X chromosome in mammals. Nat Genet. 2006;38:47–53.
Yoo MJ, Szadkowski E, Wendel JF. Homoeolog expression bias and expression level dominance in allopolyploid cotton. Heredity (Edinb). 2013;110:171–80.
Pala I, Coelho MM, Schartl M. Dosage Compensation by Gene-Copy Silencing in a Triploid Hybrid Fish. Curr Biol. 2008;18:1344–8.
Duncan CG, Grimm SA, Morgan DL, Bushel PR, Bennett BD, et al. Dosage compensation and DNA methylation landscape of the X chromosome in mouse liver. Sci Rep. 2018;8:10138.
Kelley RL, Kuroda MI. Noncoding RNA genes in dosage compensation and imprinting. Cell. 2000;103:9–12.
Wei JW, Huang K, Yang C, Kang CS. Non-coding RNAs as regulators in epigenetics (Review). Oncol Rep. 2017;37:3–9.
Tao M, Zhou Y, Li S, Zhong H, Hu H, et al. MicroRNA Alternations in the Testes Related to the Sterility of Triploid Fish. Mar Biotechnol. 2018;20:739–49.
Wang X, Werren JH, Clark AG. Allele-Specific Transcriptome and Methylome Analysis Reveals Stable Inheritance and Cis-Regulation of DNA Methylation in Nasonia. PLoS Biol. 2016;14: e1002500.
Shi X, Ng DW, Zhang C, Comai L, Ye W, et al. Cis- and trans-regulatory divergence between progenitor species determines gene-expression novelty in Arabidopsis allopolyploids. Nat Commun. 2012;3:950.
Xu P, Xu J, Liu GJ, Chen L, Zhou ZX, et al. The allotetraploid origin and asymmetrical genome evolution of the common carp Cyprinus carpio. Nat Commun. 2019;10:4625.
Ren L, Zhang X, Li J, Yan X, Gao X, et al. Diverse transcriptional patterns of homoeologous recombinant transcripts in triploid fish (Cyprinidae). Sci China Life Sci. 2021;64(9):1491–501.
Zhong H, Zhou Y, Liu S, Tao M, Long Y, et al. Elevated expressions of GH/IGF axis genes in triploid crucian carp. Gen Comp Endocrinol. 2012;178:291–300.
Yang J, Liu D, Wang X, Ji C, Cheng F, et al. The genome sequence of allopolyploid Brassica juncea and analysis of differential homoeolog gene expression influencing selection. Nat Genet. 2016;48:1225–32.
Ren L, Cui J, Wang J, Tan H, Li W, et al. Analyzing homoeolog expression provides insights into the rediploidization event in gynogenetic hybrids of Carassius auratus red var. x Cyprinus carpio. Sci Rep. 2017;7:13679.
Liu S, Luo J, Chai J, Ren L, Zhou Y, et al. Genomic incompatibilities in the diploid and tetraploid offspring of the goldfish x common carp cross. Proc Natl Acad Sci USA. 2016;113:1327–32.
Yao H, Dogra Gray A, Auger DL, Birchler JA. Genomic dosage effects on heterosis in triploid maize. Proc Natl Acad Sci USA. 2013;110:2665–9.
Wang Y, Tang H, DeBarry JD, Tan X, Li J, et al. MCScanX: a toolkit for detection and evolutionary analysis of gene synteny and collinearity. Nucleic Acids Res. 2012;40:e49–e49.
Li Q, Hermanson PJ, Springer NM. Detection of DNA Methylation by Whole-Genome Bisulfite Sequencing. Methods Mol Biol. 2018;1676:185–96.
Song Q, Zhang T, Stelly DM, Chen ZJ. Epigenomic and functional analyses reveal roles of epialleles in the loss of photoperiod sensitivity during domestication of allotetraploid cottons. Genome Biol. 2017;18:99.
Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics. 2011;27:1571–2.
Sun D, Xi Y, Rodriguez B, Park HJ, Tong P, et al. MOABS: model based analysis of bisulfite sequencing data. Genome Biol. 2014;15:R38.
Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.
Mackowiak SD. Identification of novel and known miRNAs in deep-sequencing data with miRDeep2. Curr Protoc Bioinformat Chapter. 2011. https://doi.org/10.1002/0471250953.bi1210s36.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.
Enright AJ, John B, Gaul U, Tuschl T, Sander C, et al. MicroRNA targets in Drosophila. Genome Biol. 2003;5:R1.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformat. 2014;30:2114–20.
Andrews S. FastQC: a quality control tool for high throughput sequence data. 2010
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.
Anders S, Pyl PT, Huber W. HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31:166–9.
This research was supported by the National Natural Science Foundation of China (Grant Nos. 31730098, U19A2040, 31702334), Laboratory of Lingnan Modern Agriculture Project (NT2021008), Huxiang Young Talent Project (Grant No. 2021RC3093), China Agriculture Research System of MOF and MARA (Grant No. CARS-45), 111 Project (Grant No. D20007), High-Level Talent Agglomeration Program of Hunan, China (Grant No. 2019RS1044), “Youth Talent Cultivation Program” of Hunan Normal University (XP3200202), and Natural Science Foundation of Hunan Province Grants (Grant No. 2020JJ5355).
Ethics approval and consent to participate
The animal work was approved by the academic committee in Hunan Normal University. The authors declare that animal handling complied with the relevant guidelines of the National Institutes of Health Guide for Care and Use of Laboratory Animals and the Engineering Center of Polyploidy Fish Breeding of the National Education Ministry, Hunan Normal University, Hunan, China.
The authors have declared that no competing interests exist.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
DNA methylation level of regions of TEs. A.TE distributions in each orthologous chromosomes of genomes R (2nRR) and C (2nCC). TE distributions were identified using 0.1 Mb sliding windows. B.The total methylation ratio (the combination of two homoeologs) in the two triploids and their inbred parents.C.The methylation ratio of homoeologs R vs.C in the two triploids and the total expressions in the two parents. Fig. S2. The cluster of DNA methylation of gene element among the two triploids and their inbred parents. DNA methylation in subgenomes R and C of the triploids. Black boxes indicate the DMRs, in which the similar trends of DNA methylation were in subgenomes/genomes R and C, respectively. Green box indicates the other DMRs in among these fishes. “Up in 2 k” represents 2 kb upstream of transcription start site. “Down in 2 k” represents 2 kb downstream of transcription termination site. “First exon” and “First intron” represent the first exon and intron in gene body, respectively. “Inner exon” represents all exons, except first and last exons. “Inner intron” represents all introns, except first intron. “Last exon” represents last exon in gene body. Each region was divided into twenty bins based on length. Fig. S3. GO annotation (Biological Process in level 3) of gene regulated by DNA methylation. Fig. S4. Correlation analyses between the values of differential expression (DE) and values of differential methylation (DM). A.Negative correlation between values of DE and DM (24 genes) in group of 3nR2C vs.2nCC.B.Negative correlation between values of DE and DM (5 genes) in group of 3nRC2vs.2nCC. C.Negative correlation between values of DE and DM (15 genes) in group of 3nR2C vs.2nRCC. D.Negative correlation between values of DE and DM (183 genes) in group of 3nRC2vs.2nRR. Red dot indicates the negative correlation between values of DE and DM. black dot indicates the positive correlation between values of DE and DM. Fig. S5. Differential expression (DE) analysis of miRNAs among the two triploids and their inbred parents. A.The number represents the number of up-regulated genes in the sample of the corresponding comparison. For example, the number “18” reflects that the 18 genes were higher expression in 2nRRthan in 3nRC2. B.The distribution of DE miRNAs in the five comparisons (3nR2Cvs.3nRC2, 3nR2Cvs. 2nRR, 3nRC2 vs.2nRR, 3nR2Cvs. 2nCC and 3nRC2 vs.2nCC). Fig. S6. GO annotation (Biological Process in level 3) in homoeologous recombinant genes (HRGs). Fig. S7. The distribution ofthe six DMGs in pathway of growth hormone synthesis, secretion and action (map04935). Red presents the hyper-DMGs. Green presents the HRGs in hypo-DMGs, while pink presents the other hypo-DMGs. PRKCA(K02677): classical protein kinase C alpha type; ATF4(K04374): cyclic AMP-dependent transcription factor ATF-4; CACNA1F(K04853): voltage-dependent calcium channel L type alpha-1F; PLCB(K05858): phosphatidylinositol phospholipase C, beta; CREB1(K05870): cyclic AMP-responsive element-binding protein 1; STAT5A(K11223): signal transducer and activator of transcription 5A. Fig. S8. The distribution of the nine DMGs in pathway of insulin secretion (map04911).Green presents the HRGs in hypo-DMGs, while pink presents the other hypo-DMGs. ATP1B(K01540): sodium/potassium-transporting ATPase subunit beta; PRKCA(K02677): classical protein kinase C alpha type; ATF4(K04374): cyclic AMP-dependent transcription factor ATF-4; STX1A(K04560): syntaxin 1A; CACNA1F(K04853): voltage-dependent calcium channel L type alpha-1F; KCNN2(K04943): potassium intermediate/small conductance calcium-activated channel subfamily N member 2; PLCB(K05858): phosphatidylinositol phospholipase C, beta; CREB1(K05870): cyclic AMP-responsive element-binding protein 1; SNAP25(K18211): synaptosomal-associated protein 25. Table S1. The mapping information of Methyl-seq data. Table S2. Summary of methylated and unmethylated cytosines. Table S3. The TE distribution of C. auratus red var. and C. carpio haematopterus. Table S4. The distribution of methylated and unmethylated cytosines (CpG) in the two triploids and their inbred parents. Table S5. The summary of miRNA-seq data. Table S6. The pathway distributions of hyper-DMGs and hypo-DMGs (more than 5 in each pathway).
About this article
Cite this article
Ren, L., Zhang, H., Luo, M. et al. Heterosis of growth trait regulated by DNA methylation and miRNA in allotriploid fish. Epigenetics & Chromatin 15, 19 (2022). https://doi.org/10.1186/s13072-022-00455-6