Integrated mRNA and miRNA expression profile analysis of female and male gonads in Hyriopsis cumingii

Hyriopsis cumingii is an important species for freshwater pearl cultivation in China. In terms of pearl production, males have larger pearls and better glossiness than females, but there are few reports focusing on the sex of H. cumingii. In this study, six mRNA and six microRNA (miRNA) libraries were prepared from ovaries and testes. Additionally, 28,502 differentially expressed genes (DEGs) and 32 differentially expressed miRNAs (DEMs) were identified. Compared with testis, 14,360 mRNAs and 20 miRNAs were up-regulated in ovary, 14,142 mRNAs and 12 miRNAs were down-regulated. In DEGs, the known genes related to sex determinism and/or differentiation were also identified, such as DMRT1, SOX9, SF1 for males, FOXL2 for females, and other potentially significant candidate genes. Three sex-related pathways have also been identified, which are Wnt, Notch, and TGF-beta. In 32 DEMs, the three miRNAs (miR-9-5p, miR-92, miR-184) were paid more attention, they predicted 28 target genes, which may also be candidates for sex-related miRNAs and genes. Differential miRNAs target genes analysis reveals the pathway associated with oocyte meiosis and spermatogenesis. Overall, the findings of the study provide significant insights to enhance our understanding of sex differentiation and/or sex determination mechanisms for H. cumingii.

www.nature.com/scientificreports/ of pmarg-FOXL2 and pmarg-FEM1-like for sex inversion and sex differentiation, which provided a powerful resource for the molecular mechanism of reproductive strategy in hermaphroditic mollusks 8 . These results may suggest that genes such as FOXL2, DMRT, SOX can play a role during sex regulation in many species, indicating that their functions are very conservative. MicroRNA is a single-stranded non-coding RNA molecule with a length of about 21-24 nucleotides. And it is a post-transcriptional regulator. miRNAs sequencing has advanced our understanding of miRNAs in biological sex differentiation and gonadal development, and it has been applied in a variety of aquatic species. The sex-biased miRNAs (such as miR-9, miR-21, miR-30a, miR-96, miR-200b, miR-212 and miR-7977) were obtained by a comprehensive analysis of miRNA and mRNA expression profiles in the gonads of tilapia at the early stage of sexual differentiation. Their target genes include FOXL2, AMH, STAR1, SF1 and DMRT1, which are key molecules involved in vertebrate sex differentiation 9 . Three significant miRNAs (aca-miR-30b-5p, ame-miR-263b and cfa-miR-125a) were screened from Macrobrachium nipponense. They and their predicted target genes may have a strong impact on sex differentiation/determination 10 . miRNAs were also identified and analyzed during the gonadal development of Macrobrachium rosenbergii and zebrafish 11,12 . These data will be useful for further study of miRNA mediated gonadal development and reproductive regulation mechanism. In animals, miRNAs regulate target gene expression through degradation or translation inhibition 13 . The relationship between miRNA and mRNA is not a one-to-one correspondence. A single miRNA can target multiple mRNAs, and one mRNA can also have several miRNAs binding sites 14 . During gonadal development, miRNAs are expressed differently in male and female gonads. They play a role through germ cells and gonadal cells to regulate critical proteins needed for gonadal development 15 . Studies have shown that miR-202-5p/3p can be expressed in a sexual dimorphism pattern as the primordial XY gonad differentiates into a testis and play an early role in mouse testis development 16 . The results of knockdown or overexpression of miR-124 in mouse showed that miR-124 could induce the repression of both SOX9 translation and transcription in ovarian cells 17 . Therefore, understanding the regulatory relationship between miRNA and sex-determining genes can help us have a better comprehending of the sex-determining mechanism. Bioinformatics analysis of miRNA and mRNA expression profile can help us improve the reliability of prediction and understand the molecular mechanism of post-transcriptional regulation. Nowadays, there is no information about the interaction network and regulation of mRNAs and miRNAs in the gonads of H. cumingii. Therefore, our goal is to identify the essential genes and miRNAs that regulate sex determination/differentiation of H. cumingii by mRNA sequencing and miRNA sequencing. This study provides basic information for understanding sex differentiation and sex determining mRNAs and miRNAs of H. cumingii.

Results
mRNA sequencing and assembly. The samples were sequenced from the gonads of six H. cumingii, including three ovaries (F1, F2, and F3) and three testes (M1, M2, and M3). A total of 308,008,006 raw reads were produced by Illumina HiSeq X Ten sequencing. And 302,865,864 clean reads filtered by Trimomatic. The average percentage of Q30 and the G + C percentage of the six cDNA libraries were 94.75% and 41.48% (Table 1) Screening for DEGs and functional annotation. Three female gonadal tissues (F1 ~ F3) were Group_F, and three male gonadal tissues (M1 ~ M3) were Group_M. Compared with Group_M, there were 28,502 DEGs in Group_F, including 14,360 genes up-regulated and 14,142 genes down-regulated. GO annotations consist of three categories: biological processes, cellular components, and molecular functions. DEGs were enriched to 9295 GO terms. Ten terms with P value less than 0.05 and the largest number of DEGs were selected in each of the three categories ( Supplementary Fig. 1). KEGG enrichment of DEGs showed that there were 305 pathways, remove the pathways in which the number of DEGs is less than three, sort them according to the − log 10 P value, and draw a bubble diagram with the first 20 pathways (Supplementary Fig. 2).
According to sex determining genes of other species 18,19 and GO analysis, we screened 12 sex-related genes in H. cumingii transcriptomes (Table 2), including DMRT, TRA , SOX, FOXL genes, etc. GO analysis is to find candidate genes by searching for "sex" keywords in GO terms. The search results include nine GO terms, namely Table 1. Summary information and analysis of mRNA sequences in H. cumingii. Raw reads: the statistics of the raw reads; Clean reads: remove the contaminated and low quality read; Q30: The ratio of bases with a value greater than 30 to the total bases in raw bases; GC: The ratio of the total number of G and C to the total number of bases in clean bases.   Table 2. In the KEGG enrichment results, Wnt (ko04310), Notch (ko04330), and TGF-beta (ko04350) were confirmed to be related to sex determination 20 . There are 49 DEGs in these three signaling pathways, which may be related to sex regulation of H. cumingii (Supplementary Table 1). Through these methods, we obtained a total of 67 sex-related candidate genes (after removing duplicates). 14 differentially expressed genes were selected for quantitative real-time PCR (qRT-PCR) detection ( Supplementary Fig. 3).
miRNAs sequencing and assembly. We constructed six cDNA libraries of small RNAs from F1, F2, F3, M1, M2, and M3. After filtering the Rfam database, transcript sequence, and Repbase database, the reads are compared with miRBase to identify and annotate the known miRNA (  Table 2). Among the known miRNAs, three were ovary-biased, and ten were testis-biased. Unsupervised hierarchical clustering groups the samples in two major clusters (Fig. 1). The number of target genes of 32 DEMs are 911. The functional terms cytoplasm, ATP binding, and nucleoplasm, had a higher number of transcripts than other terms ( Supplementary Fig. 4). KEGG analysis showed that the target genes were enriched in 184 signaling pathways. The largest number of target genes is found in cell cycle (Supplementary Fig. 5). To study the expression pattern of these DEMs in testis and ovary, we randomly selected three ovary-biased miRNAs (miR-184, miR-281-3p, miR-315), and five testis-biased miRNAs (miR-153, miR-29, miR-7, miR-9-5p, miR-92) to verify their expression profiles by qRT-PCR (Fig. 2). All of these miRNAs expression patterns are consistent with the results of miRNAs sequencing, indicating that the sequencing results are reliable.

Integrated analysis of DEMs and DEGs.
According to the results of transcriptome sequencing and miR-NAs sequencing, we found three key miRNAs (miR-9-5p, miR-92, miR-184) and their target genes (Table 4). These genes may be related to the sex regulation of H. cumingii. Some of them have been reported to relates to  www.nature.com/scientificreports/  www.nature.com/scientificreports/ sex. For example, REX1 is associated with X-chromosome inactivation in mice 21 . SMAD4 is essential factors that regulate the female fate of human germ cells 22 . MCM6 was identified and considered as candidate genes that may be involved in sex differentiation regulation in cucumber 23 . CDH5 is necessary for spermatogenesis and chromatin concentration in mice 24 . Besides, 361 differentially-expressed target genes were obtained from the intersection of 28,502 DEGs and 911 target genes (Fig. 3a). GO analysis showed that there were more transcripts in the cytoplasm, ATP binding, calcium ion binding, endoplasmic reticulum membrane, and microtubule (Fig. 3c).

Discussion
In order to explore the mRNAs and miRNAs involved in sex differentiation and/or sex determination of H. cumingii, RNA-seq and miRNA-seq analyses were performed on ovaries and testes. To the best of our knowledge, this is the first report of miRNA and mRNA profiling of ovary and testis in H. cumingii. We analyzed miRNAs and mRNAs to provide a basis for screening candidate miRNAs and sex-related genes. This result provides an essential insight into the mechanism of sex differentiation/determination in H. cumingii.
Since there are few studies focusing on the sex of bivalve shellfish, to find candidate genes related to sex determination, we refer to other species genes that have been identified, which has important reference value for our research. After screening the transcriptome results, we obtained seven DEGs as SOX9, DMRT1, FOXL2, SF1, TRA2, BMP2 and ZGLP1. In general, SOX9 is a direct target of SRY gene in mammals and is necessary for normal testicular development 25 . The conservative expression of SOX9 in chicken and mouse indicates its primary role in the testicular determination of vertebrates 26 . SOX9 was also identified in bivalve shellfish, such as N. subnodosus, C. gigas, and Hyriopsis schlegelii 4,27,28 . In our results, miranda software predicted that there was an interaction site between miR-193 and SOX9. Also, SOX9 was highly expressed in testis ( Supplementary Fig. 3). These results suggest the potential role of SOX9 and miR-193 in the sex of H. cumingii. Unfortunately, due to the limitation of miRNA sequencing information, other genes were not predicted by miRNAs except SOX9. Apart from that, DMRT (double sex/male-normal-3-related transcription factor) gene is a major transcription factor controlling sex determination and differentiation 29 . DMRT1 can not only control testicular differentiation and meiosis of male germ cells in mice 30,31 . But also be necessary for male sexual development in zebrafish 32 . In mollusks, the homologous gene of DMRT has been identified in many species, such as C. gigas, P. margaritifera, and P. yessoensis [33][34][35] . In H. cumingii, the qRT-PCR results showed that DMRT1 gene was a testis specific gene, and its expression level in the ovary was very low, and the expression level in mRNA sequencing results was not detected (Supplementary Fig. 3). It was suggested that DMRT1 was crucial for the male sexual development of H. cumingii. FOXL2, a member of the forkhead box (Fox) domain transcription factor family, plays a key role in ovarian differentiation and oogenesis in vertebrates 36 . FOXL2 can affect ovarian development and sex determination of mice 37 . In C. farreri and P. yessoensis, FOXL2 expression was a sexually dimorphic pattern, and the expression level in ovary was significantly higher than that in testis 38,39 . In our results, FOXL2 was highly expressed in the ovary, and the expression abundance in the testis was very low. The results of mRNA sequencing also showed that FOXL2 was not expressed in testis (Supplementary Fig. 3). It is suggested that FOXL2 is crucial for the female sexual development of H. cumingii. SF1 (Steroidogenic factor 1) is an important transcription factor involved in steroidogenesis, reproduction, and sexual differentiation 40 . The TRA2 gene, encoding a protein (TRA-2) which directs sex-specifically alternative splicing, has been proved to play important roles in sex differentiation and sex development of Drosophila melanogaster 41 . BMP2 (Bone morphogenetic protein 2) has been confirmed to be related to the fetal ovarian development and gonadal somatic cell differentiation of mice 42,43 . ZGLP1 (GATAtype zinc finger protein 1) is essential for the oogenic program and meiotic entry in mice 44 . In addition to the above genes, GO analysis also enriched some extra sex determination/differentiation genes, they are CYP17A1, DACH2, DPN, SMC5, and SRD5A1. CYP17A1 has the activity of 17α-hydroxylase and 17, 20-lyase, which is www.nature.com/scientificreports/ involved in the steroidogenic pathway that produces androgens and estrogens, it is crucial for male determination in amphibians 45 . DACH is differentially expressed in the male and female genital discs of drosophila, and plays sex-specific roles in the developmenting genitalia 46 . And DPN is also a primary part of sex determining signal of drosophila 47 . SMC5 was confirmed to be associated with spermatogenesis of mice 48 . In Crassostrea hongkongensis, SRD5A1 was expressed at a higher level in female gonads than in other tissues 18 . There is no doubt that the above role in gender regulation is vital. In the sequencing results, FOXL2, BMP2, ZGLP1, CYP17A1, DACH2, DPN, SRD5A1 were highly expressed in the ovary, and DMRT1, SF1, SMC5, SOX9 were highly expressed in the testis. This suggests that they may play a role in the formation of female or male. The results of random verification of qRT-PCR showed the accuracy of sequencing ( Supplementary Fig. 3), but their position in sex regulation needs more research to explore. In KEGG results, the function of Wnt, Notch, and TGF-beta signaling pathways in sex determination and gonadal development in mammals have been relatively straightforward, including BMPs (bone morphogenetic proteins), HH (hedgehog), Hippo, nuclear receptors 20,49 . In addition to www.nature.com/scientificreports/ non-mammals, it has also been shown that the classical Wnt signaling pathway regulates gonadal differentiation in zebrafish 50 . The Notch signaling pathway can regulate stem cell proliferation and sex determination in the germline of Caenorhabditis elegans 51 . The role of AMH (anti-Müllerian hormone) in testicular differentiation of Eriocheir sinensis highlights the importance of TGF-beta pathway in sex determination of reptiles 52 . All three pathways exist in our results, indicating that they may be related to the sex determination/differentiation of H. cumingii. At the same time, in these three signal pathways, some gender-related genes have also been studied. WNT2, a wnt family ligand, has been demonstrated to be involved in the induction of sex-specific differences in drosophila 53 . HES1 isan important transcription factor in Notch signaling pathway, which is related to the mouse testis cell fate determination step from primordial stage to the differentiated stage in adulthood 54 . The loss of CREBBP (CREB binding protein, Notch and TGF-beta signaling pathway member) and P300 (its related paralogue) disrupts histone acetylation of mouse SRY promoter and causes XY gonadal sex reversal 55 . CUL-1 (CULLIN-1, TGF-beta signaling pathway member), SKR-1, and SEL-10 constitute a SCF E3 ligase complex that plays an critical role in modulating sex-determination of C. elegans 56 . In this study, 32 miRNAs were identified from six miRNA libraries, among which 19 were novel miRNAs. Since novel miRNAs have not been studied in other species, we focus on 13 known miRNAs. The DEMs in our gonadal miRNA-seq results are similar to those of other species. For example, miR-7, miR-9, miR-29, miR-184 has been detected in miRNA-seq identification of other species [57][58][59] . In these miRNAs, miR-7 is a crucial factor in regulating the hypothalamus-pituitary-ovary axis of mice 60 . Moreover, miR-7a2 can regulate sexual maturation and reproductive function 61 . miR-9 is involved in spermatogenesis of spermatogonia during the natural sex change of Monopterus albus 62 . Furthermore, it can also regulate the critical genes in the ovarian development pathway of Mud Crab Scylla paramamosain 63 . miR-29 is an essential regulatory factor in the process of male meiosis of mice 64 . miR-184 was highly expressed in mouse testis, and down-regulated the expression of NCOR2 during spermatogenesis 65 . In drosophila, miR-184 can affect oogenesis, early embryonic development, and ovarian development 66,67 . DROSHA is essential for the coordinated development of somatic and germ cell precursors in drosophila larvae. And miR-8 and miR-184 are regulated by DROSHA 67 . miR-92 may play a role in early chicken gonadogenesis 68 . In Japanese flounder, miR-92 regulates early development and metamorphosis of Paralichthys olivaceus 69 . We list three miRNAs and their target genes, and there are 28 DEGs (Table 4). Based on the potential function of these miRNAs in sex, it is speculated that these DEGs may also be sex-related candidate genes of H. cumingii. However, the specific functionality needs further experimental verification. Also, qRT-PCR results showed that miR-7, miR-9-5p, miR-29, miR-92, and miR-153 were highly expressed in the testis and miR-184, miR-218, miR-315 were highly expressed in the ovary (Fig. 2). This is consistent with our sequencing results, which further indicates that these DEMs may play a specific role in the sex regulation of H. cumingii. 361 differentially-expressed target genes were obtained from the intersection analysis of 28,502 DEGs genes and 911 target genes (Fig. 3a). These genes are differentially expressed in ovaries and testes, and their associated miRNAs are also differentially expressed between the two genders. There are two terms of spermatogenesis and oocyte meiosis in GO and KEGG analysis (Fig. 3b,c). They have seven differentially expressed genes (IFT81, LRGUK, SYNE1, TDRD1, PPP1C, PKMYT, and CYCE). In mouse spermatogenesis, IFT81 plays an essential role in regulating the assembly and elongation of sperm flagellum 70 . And LRGUK-1 is required for basal body and manchette function during spermatogenesis and male fertility 71 . PPP1C is needed in the meiosis of mouse oocytes 72 . In P. olivaceus, PoTDRD1 is a germ line specific and sex dimorphic factor, which may play a role in the development of its reproductive system and gametogenesis 73 . In Xenopus, CPEB and miR-15/16 co-regulate CYCE1 and play a role in oocyte maturation 74 . These genes and miRNAs may be closely related to germ cell maturation and gametogenesis of H. cumingii.

Conclusion
In this study, the transcriptome and miRNA analysis were applied to explore candidate genes/miRNAs for sex determination/differentiation in H. cumingii. mRNA and miRNA expression profiles provided a rich list of genes and miRNAs expressed in testes and ovaries. In total, 95 DEGs, three DEMs, and three signaling pathways were screened, which may be closely relative with sex determination/differentiation in H. cumingii. In the intersection analysis of DEGs and DEMs target genes, two signal pathways related to gametogenesis were obtained. They include seven genes and seven miRNAs. These sex-related genes, miRNAs, and signaling pathways provided not only the basis foundation for the study of the sexual regulation mechanism in H. cumingii, but also a fundamental reference for better understanding the sex determination and/or differentiation of bivalves.

Materials and methods
Ethics statement. The samples were collected from Jinhua farm in Zhejiang province. All experimental protocols were approved by the Institutional Animal Care and Use Committee (IACUC) of Shanghai Ocean University, Shanghai, China.
Sample collection. Specimens used in current study were obtained WeiMing aquaculture farms (ZheJiang, China). Three sexually mature females and three males H. cumingii weighing 320-338 g were collected. In general, the mussels over 2 years old are considered to be sexually mature. The samples we used in this study were 3 years old. Because the eggs of H. cumingii are very large and the sperm is very small. We took a small amount of liquid from the gonads and observed them under a microscope to judge their sex. These mussels are propagated from the same batch of parents and cultured in the same culture pond. The water temperature in the pond ranges from 12 to 32 °C throughout the year, and the pH is about 7.5. All the samples were taken back to the laboratory of Shanghai Ocean University. Keep the mussels in a tank for 3 days. The water temperature was controlled around 26 °C. The gonadal tissues of three female mussels (F1, F2, F3) were Group_F, and those of three male www.nature.com/scientificreports/ mussels (M1, M2, M3) were Group_M. Gonadal samples were collected with high temperature sterilized scissors and tweezers, quickly placed in liquid nitrogen, and then transferred to the − 80 °C refrigerator. cDNA library and small RNA preparation and sequencing. Six sequencing libraries were constructed, three libraries (F1, F2, F3) from female groups and other three (M1, M2, M3) from male groups. The samples were prepared by using a mirVana miRNA Isolation Kit (Invitrogen, USA). The libraries were constructed by using TruSeq Stranded mRNA Library Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. After passing the quality test of Agilent 2100 Bioanalyzer, the constructed library was sequenced by Illumina HiSeq X Ten sequencer to produce the double-terminal data of 150 bp. The miRNA library was prepared by using mirVana miRNA Isolation Kit (Invitrogen, USA). The main experimental steps were the ligation of 3′ and 5′ connectors. The small RNAs equipped with connectors were reversed transcript and amplified by PCR. Then 147 nt and 157 nt bands were recovered by RNA gel electrophoresis. The concentration of RNA was detected by Agilent 2100 Bioanalyzer. The platform of small RNA sequencing was Illumina HiSeq 4000. Transcriptome and small RNA sequencing were both completed by OE biotech company (Shanghai).

Assembly and annotation.
A large number of mRNA-seq raw reads were obtained from male and female gonadal tissues. Remove adaptor, low quality reads and containing ploy-N reads from the raw reads. The clean reads were assembled into expressed sequence tag clusters (contigs). The transcripts were assembled with Trinity (vesion 2.4) 75 . The longest transcript was selected as a unigene based the similarity and length of a sequence analysis. The transcriptome analysis was conducted by OE biotech company (Shanghai, China). The transcript annotations give functional annotations for NR (ftp://ftp.ncbi.nih.gov/blast /db), SwissProt (http://www.unipr ot.org/downl oads), Clusters of orthologous groups for eukaryotic complete genomes (KOG) (ftp://ftp.ncbi.nih. gov/pub/COG/KOG/kyva), Pfam, GO (Gene Ontology) classification, and KEGG (Kyoto Encyclopedia of Genes and Genomes). The Pfam database was compared with the protein family model by HMMER 76 software. Other databases were annotated by Diamond 77 software with a threshold E-value of 10 −5 . Base on the SwissProt annotation, GO classification was performed by the mapping relation between SwissProt and GO term. The unigenes were mapped to the KEGG 78,79 database to annotate their potential metabolic pathways.
The small RNA sequencing reads were converted into sequence data (also called raw data/reads) by base calling. After filtering out the reads without 3′ adapter and insert tag, reading shorter than 15 nt and longer than 41nt, the clean reads were obtained. The clean reads sequences were compared with Rfam v.10.1 (http://www. sange r.ac.uk/softw are/Rfam) by blastn. rRNA, scRNA, Cis-reg, snRNA, tRNA and other sequences were annotated. And the transcripts and repetitive sequences were filtered. The filtered sequence was used to identify the known miRNAs by miRBase v.21 database (http://www.mirba se.org/). Unannotated small RNAs were analyzed by mirDeep2 80 . The targets of differentially expressed miRNAs were predicted by using software miRanda 81 , with the parameter as follows: S ≥ 150 ΔG ≤ − 30 kcal/mol and demand strict 5′ seed pairing.
DEGs and DEMs analysis, GO and KEGG enrichment. Bowtie2 and eXpress were applied to analyze FPKM (fragments per kilobaseper million mapped reads) of each unigene. The counts number of each sample unigene was standardized by using DESeq software. The fold change was calculated. NB (negative binomial distribution test) was used to test the difference significance of reads number. P value adjusted with a falsediscovery rate (FDR) correction for multiple testing by Benjamini-Hochberg method 82 . Finally, the DEGs were screened according to the |fold change|> 2 and P < 0.01. After getting the DEGs, GO and KEGG enrichment analysis were carried out to determine the biological functions or pathways. The number of differential mRNAs included in each GO term was counted. The significance of differential gene enrichment in each GO term was calculated by hypergeometric distribution test. The result of the calculation will return a P value of enrichment significance. And the small P indicates that the differential unigene is enriched in the GO term. Similarly, KEGG database was used to for analysis DEGs. The hypergeometric distribution test method was used to calculate the significance of differential unigene enrichment in each Pathway. At the same time, the unsupervised hierarchical clustering of DEGs is carried out. And the expression pattern of DEGs among different samples is shown in the form of heat map. miRNA expression is calculated using TPM (transcript per million), TPM = the number of read compared to each miRNA / the total ratio of samples to the number of read × 10 6 . The miRNAs with P value < 0.05 and fold change more than 2 times were selected as differentially expressed miRNAs. Using unsupervised hierarchical clustering analysis, a heat map is constructed based on differentially expressed miRNAs. The heat map (Fig. 1) is drawn by the pheatmap package in R 83 . GO enrichment and KEGG pathway enrichment analysis of different expressed miRNA-target-Gene were respectively performed using R based on the hypergeometric distribution. Figure 3a is drawn by VennDiagram package in R 84 . Figure 3b,c are drawn by ggplot2 R package 85 . qRT-PCR of mRNA and miRNAs. We selected 12 genes and 8 miRNAs to verify the sequencing results by qRT-PCR. The primer sequences of mRNA and miRNA are shown in supplementary Table 3 and supplementary  Table 4, respectively. For the verification of differentially expressed genes, the kit for RNA reverse transcription is PrimeScript RT reagent Kit (TaKaRa, Japan). The total qRT-PCR reaction of mRNA contained 2 × TB Green Pre- www.nature.com/scientificreports/ mix Ex Taq II (TaKaRa, Japan) 10 μL, forward and reverse gene-specific primer (10 mM) 0.8 μL, cDNA template 1.6 μL, ddH 2 O 6.8 μL. The condition of the reaction is 95 °C (pre-denaturation) for 15 min, 95 °C (degeneration) for 10 s, and 60 °C (annealing) for 30 s (40 cycles), followed by dissociation curve analysis at 95 °C for 15 s, 60 °C for 1 min and 95 °C for 15 s. EF1α was used as the internal reference gene. The reverse transcription of miRNAs was performed using the Mir-X miRNA First-Strand Synthesis kit (TaKaRa, Japan). The total qRT-PCR reaction of miRNAs contained ddH 2 O 9 μL, TB green advantage premix 12.5 μL, Rox dye 0.5 μL, miRNA specific primer (or U6 forward primer) 0.5 μL, MRQ 3′primer (or U6 reverse primer) 0.5 μL, and cDNA template 2.0 μL. The reaction procedure was set according to the following procedure: 95 °C for 10 s; 95 °C for 5 s, 60 °C for 20 s, 39 cycles; then the PCR temperature was increased from 60 to 95 °C to generate the dissolution curve. U6 snRNA was used as internal reference gene. All qRT-PCR were performed using Bio-Rad CFX-96 (Bio-Rad, USA). The relative expression of mRNAs and miRNAs was calculated by 2 −ΔΔCT method. The difference was calculated by SPSS12.0. Statistically significant differences were examined by paired t test. A value of P < 0.05 was considered to statistically significant.

Data availability
The datasets generated and analyzed during this are available from the corresponding author on reasonable request. www.nature.com/scientificreports/