Morphological changes of strobili during the flower development
To understand the morphogenesis of female and male strobili of G. biloba and their important functions during pollination, the changes in macroscopic and ultrastructure during their development was observed in this study. As shown in Fig. 1, male flower buds are borne at the apex of the short shoots had a dormancy stage of about 6 months. After the dormancy stage, male buds gradually unfolded and mature (Fig. 1ae). Each male cone contained 6080 microsporophylls, and each microsporophyll is composed of a sterile extension and two elliptical microsporangia (Fig. 1fg). Female flower buds are mainly surrounded by thick bud scales until mid-late March, when bud scales began to unfold afterwards. By mid-April, each ovule consisted of a slender stem, an ovule ring, an integument and a nucellus during pollination (Fig. 1hk). The results of scanning electron microscopy (SEM) showed that ovules differentiated into many structures, including the micropyle, micropylar canal and pollen chamber before pollination (Fig. 1im).
Comparison of phytohormone contents of female and male strobili of G. biloba
Zeatin (ZT), GA3 and IAA were detected at different developmental stages of the female and malestrobili of G. biloba (Fig. 2). It was worth mentioned that the change trend of hormone contents was consistent in two developmental stages of male and female strobili. In detail, the contents of ZT and IAA in male strobili were significantly higher than those in female strobili, while the contents of GA3 in female strobili were significantly higher than those in male strobili (Fig. 2).
Overview of small RNA deep-sequencing data
For identifying miRNAs involved in the sexual differentiation and development of flower in G. biloba, four small RNA libraries of MB, FB, MS and OS were constructed. A total of approximately 79.489 million raw reads were produced. After removing contaminated reads, we obtained clean reads ≥10.78 M from each sample. A total of 3,293 miRNAs were obtained from all samples, of which 1,085 were conserved miRNAs and 2,208 were novel miRNA (Additional file 1: Table S1). All miRNAs sequences were in the range of 1825 nt, mainly concentrated in 2024 nt, of which 21 nt was the most abundant type in the small ribonucleic acid library (Fig. 3A). The above-mentioned results are consistent with the typical length of mature miRNA in other plant species [30]. The first base of the 5’ end of the mature miRNA has a strong bias toward U (Fig. 3B). In all base sites, the proportion of bases A and U exceeds 50% except for sites 18 and 19 (Fig. 3C). These statistical results are basically consistent with the characteristics of plant miRNAs, indicating the accuracy of high-throughput sequencing results.
Differential expression of miRNAs
To study the mechanism of sexual differentiation in G. biloba strobili, a hierarchical cluster analysis of the differentially expressed miRNAs was performed for grouping miRNAs based on expression level (Fig. 4A). Among 239 differentially expressed miRNAs between MB _vs_ FB, 103 were up-regulated and 136 were down-regulated in FB (Fig. 4B). A total of 247 miRNAs were differentially expressed between MS _vs_ OS, of which 142 miRNAs were up-regulated and 105 miRNAs were down-regulated in FB (Additional file 2: Table S2). For analyzing the miRNAs involved in the regulation of sexual differentiation in G. biloba strobili, this study mainly analyzed the differentially expressed miRNAs in the same stage. There were 61 miRNAs continuously down-regulated in male strobili (Fig. 4C) and 36 continuously up-regulated in female strobili (Fig. 4D).
Functional annotation of target genes of differentially expressed miRNA
To investigate the regulatory function of miRNA with respect to target genes, functional annotation was performed on target genes corresponding to differentially expressed miRNA. A total of 3,043 target genes out of 6,032 genes were annotated to 45 GO categories composed of 13 (2,450 target genes, 80.51%), 11 (1,308 target genes, 42.98%) and 20 (2,129 target genes, 69.96%) categories, which belonged to the molecular function, cellular components, and biological processes, respectively (Fig. 5A). For extending the understanding of the biological functions at the protein level, a total of 2,395 target genes were annotated to 25 protein families with classified COG functions (Fig. 5B). Among them, the “general function prediction only” accounted for the largest proportion (627 target genes), followed by “replication, recombination, and repair” (492 target genes). A total of 1,251 target genes were successfully annotated to 97 KEGG pathways, and the top 50 metabolic pathways with a higher number of annotated genes were shown in Fig. 5C. In particularly, 60, 59 and 48 target genes were involved in plant hormone signal transduction (ko04075), biosynthesis of amino acids (ko01230) and phenylpropanoid biosynthesis (ko00940), respectively.
Screening of differentially expressed genes (DEGs)
A total of 4,346 DEGs was screened between MB and FB. Of which, 1,971 and 2,375 DEGs were up-regulated and down-regulated in FB _vs_ MB, respectively (Fig. 5D). By comparison MS and OS, 7,087 DEGs were obtained with 3,190 up-regulated and 3,897 down-regulated DEGs in OS. These DEGs may play role in phenotype differentiation of male and female strobili in G. biloba.
Regulatory analysis between miRNAs and mRNA
Using differentially expressed miRNAs as screening criteria, differentially expressed mRNAs regulated by differentially expressed miRNAs were obtained. The results showed 424 and 493 pairs of miRNA-mRNA combinations with up-regulated and down-regulated miRNA in MB, respectively (Table 1). Among them, 445 pairs had negative correlation. A total of 840 miRNA-mRNA regulatory pairs with up-regulated miRNAs were found in OS, and 776 miRNA down-regulated pairs were observed in OS (Table 1). A total of 846 pairs with negative regulatory relationship was screened in MS _vs_ OS (Table 1).
Table 1 Number of pairs between differentially expressed miRNAs with their target genes
DEG_Set
|
Total
|
Up
|
Down
|
MB _vs_ FB
|
4,224
|
424
|
493
|
MS _vs_ OS
|
5,050
|
840
|
776
|
Note: Up, down indicates how the miRNA-regulated target genes were regulated in this group comparison.
Identification of differentially expressed miRNAs and target genes (DEGs) involved in regulating sexual differentiation in G. biloba
Many factors play a crucial role in the process of sexual determination. Phytohormones has been proved to be crucial for plant sexual differentiation [31]. Based on the annotation information of target genes and the negative regulatory relationship between miRNAs and mRNAs, we screened candidate miRNA-regulated target genes related to the differential expression of sexual determination in G. biloba. As shown in Table 2, five miRNA-mRNA negative correlation pairs were found in the flower development pathway, two pairs of miRNA-mRNA in MB_vs_ FB corresponded to two target genes (GIGANTEA, POX53),, and three pairs of miRNA-mRNA in MS _vs_ OS corresponded to two target genes (ATR and AP2).. In MB_vs_ FB, novel_miR_1881 argets the GIGANTEA and up-regulated in FB, novel_miR_2912 argets the POX53 and down-regulated in FB. In MS _vs_ OS, novel_miR_1858 targets the ATR gene, novel_miR_1944 and novel_miR_2329 simultaneously target the AP2 gene, and were down-regulated in OS.
Table 2 MiRNAs and target genes involved in regulating the process of G. biloba strobilus development and sex differentiation
Metabolic pathway /Development
|
MB _vs_ FB
|
MS _vs_ OS
|
miRNAID
|
miR_regulate
|
mRNATar
|
Tar_regulate
|
miRNAID
|
miR_regulate
|
mRNATar
|
Tar_regulate
|
Flower development pathway
|
novel_miR_1881
|
up
|
GIGANTEA
|
down
|
novel_miR_1858
|
down
|
ATR
|
up
|
novel_miR_2912
|
down
|
POX53
|
up
|
novel_miR_2329
|
down
|
AP2
|
up
|
|
|
|
|
novel_miR_1944
|
down
|
AP2
|
up
|
GA metabolic pathway
|
miR159.1
|
up
|
GAMYB1
|
down
|
miR159.1
|
up
|
GAMYB1
|
down
|
miR159.1
|
up
|
GAMYB2
|
down
|
miR159.1
|
up
|
SRG1
|
down
|
miR159.2
|
up
|
GAMYB3
|
down
|
miR159.1
|
up
|
GAMYB4
|
down
|
miR159.2
|
up
|
GAMYB4
|
down
|
miR159.2
|
up
|
GAMYB4
|
down
|
miR159.2
|
up
|
GAMYB1
|
down
|
miR159.2
|
up
|
GAMYB1
|
down
|
miR159.3
|
up
|
GAMYB4
|
down
|
miR159.3
|
up
|
GAMYB4
|
down
|
miR159.3
|
up
|
GAMYB1
|
down
|
miR159.3
|
up
|
GAMYB1
|
down
|
miR159.4
|
up
|
GAMYB4
|
down
|
miR159.4
|
up
|
GAMYB4
|
down
|
miR159.4
|
up
|
GAMYB1
|
down
|
miR159.4
|
up
|
GAMYB1
|
down
|
miR858
|
up
|
GAMYB1
|
down
|
miR858
|
up
|
GAMYB1
|
down
|
miR858
|
up
|
GAMYB3
|
down
|
miR858
|
up
|
GAMYB4
|
down
|
miR858
|
up
|
GAMYB4
|
down
|
novel_miR_2221
|
up
|
GALMADRAFT
|
down
|
miR319.1
|
up
|
GAMYB3
|
down
|
novel_miR_3455
|
up
|
GALMADRAFT
|
down
|
miR319.2
|
up
|
GAMYB3
|
down
|
novel_miR_1944
|
down
|
PME13
|
up
|
miR319.3
|
up
|
GAMYB3
|
down
|
novel_miR_2329
|
down
|
PME13
|
up
|
novel_miR_1872
|
up
|
ARHGAP7
|
down
|
|
|
|
|
ETH metabolic pathway
|
miR2950.1
|
up
|
ERF2
|
down
|
|
|
|
|
miR2950.2
|
up
|
ERF2
|
down
|
|
|
|
|
IAA metabolic pathway
|
|
|
|
|
miR160.1
|
up
|
ARF18.1
|
down
|
|
|
|
|
miR160.1
|
up
|
ARF18.2
|
down
|
|
|
|
|
miR160.2
|
up
|
ARF18.1
|
down
|
|
|
|
|
miR160.2
|
up
|
ARF18.2
|
down
|
|
|
|
|
miR160.3
|
up
|
ARF18.1
|
down
|
|
|
|
|
miR160.3
|
up
|
ARF18.2
|
down
|
|
|
|
|
miR160.4
|
up
|
ARF18.1
|
down
|
|
|
|
|
miR160.4
|
up
|
ARF18.2
|
down
|
|
|
|
|
miR160.5
|
up
|
ARF18.1
|
down
|
|
|
|
|
miR160.5
|
up
|
ARF18.2
|
down
|
|
|
|
|
miR160.6
|
up
|
ARF18.1
|
down
|
|
|
|
|
miR160.6
|
up
|
ARF18.2
|
down
|
CTK metabolic pathway
|
|
|
|
|
miR2950.3
|
up
|
AHP5
|
down
|
A total of 23 pairs of negative regulatory relationships were involved in GA metabolic pathway, 16 pairs and 5 target genes (GAMYB1, GAMYB2, GAMYB3, GAMYB4, ARHGAP7) in the MB _vs_ FB; 15 pairs and 5 target genes (GAMYB1, SRG1, GAMYB4, GALMADRAFT, PME13) in the MS _vs_ OS. In addition, eight conserved miRNAs (miR159.1, miR159.2, miR159.3, miR159.4, miR858, miR319.1, miR319.2, and miR319.3) and novel_miR_1872 were up-regulated expressed level of observed in the FB; two novel miRNAs (novel_miR_2221 and novel_miR_3455) were up-regulated expressed level, and novel_miR_1944 and novel_miR_2329 were down-regulated in OS (Table 2). The corresponding target genes were transcription factor GAMYB (GAMYB1, GAMYB2, GAMYB3, GAMYB4) and other genes involved in plant GA metabolic pathway (ARHGAP7, GALMADRAFT),, which were down-regulated; however, the PME13 was up-regulated in OS.
For the ETH metabolic pathway, two miRNA-mRNA pairs were observed to negatively correlate in MB _vs_ FB, corresponding to one target gene (ERF2),, and down-regulated expression in FB.
In the MS _vs_ OS, twelve pairs annotated to the metabolic pathways of IAA and two target genes (ARF18.1, ARF18.2).. MiRNA160.1, miRNA160.2, miRNA160.3, miRNA160.4, miRNA160.5, miRNA160.6 were up-regulated in OS. The corresponding target genes were auxin response factor 18 (ARF18.1, ARF18.2),, which were down-regulated.
According to the analysis of transcriptome data, 235 DEGs were annotated to the plant hormone signal transduction pathway, of which 81 DEGs were in the MB _vs_ FB and MS _vs_ OS. The expression patterns of genes annotated to IAA, GA, ETH, CTK and the ABA signal transduction pathway were analyzed (Fig. 6) to find out that there were 52 DEGs, of which 29 DEGs were in the IAA signal transduction pathway (Fig. 6A), 7 DEGs in the GA signal transduction pathway (Fig. 6B), 4 DEGs in the ETH signal transduction pathway (Fig. 6C), 8 DEGs in the ABA signal transduction pathway (Fig. 6D), and 5 DEGs in the CTK signal transduction pathway (Fig. 6E). These DEGs involved in hormone signal transduction may participate in sexual determination of G. biloba. Combined with the negative correlation regulation between miRNA and mRNA, only 1 of the 52 DEGs had a negative correlation regulation, namely, AHP5 involved in the CTK signal transduction pathway. In OS, the gene was down-regulated while the corresponding miRNA (miR2950.3) was up-regulated (Table 2). This miRNA may play an important regulatory role in flower differentiation and sexual differentiation of G. biloba.
Confirmation via degradome approaches
A total of 2.283 million clean data were obtained by degradome sequencing from G. biloba strobili. Based on conserved miRNAs library, the prediction of miRNAs in the small RNA sequencing database, gene transcript sequence information files of corresponding species, and degradation site detection was performed by Cleaveland software [32]. The analysis revealed 155 degraded target genes and 156 target gene degradation sites. According to the integration analysis of degradome sequencing and the previous analysis results (Additional file 3: Table S3), fifteen miRNA-mRNA regulatory relationships were verified from the degradome sequencing, including three target genes and nine miRNAs (Fig. 7). MiR159.2, miR159.3 and miR159.4 jointly regulated GAMYB1 involved in the GA metabolic pathway. MiR160.1, miR160.2, miR160.3, miR160.4, miR160.5 and miR160.6 jointly regulated the ARF18.1 and ARF18.2 involved in the metabolic pathway of IAA (Fig. 8). There was only one miRNA degradation site on the three target genes, and the degradation site of GAMYB1 was 994, the site of ARF18.1 was 1439, and the site of ARF18.2 was 1508. These results verified that multiple miRNAs could simultaneously regulate the same target gene.
Validation of miRNA and its target gene by qRT-PCR
For verifying the identified expression patterns of miRNAs and corresponding target genes involved in the process of sexual differentiationof G. biloba strobili, 12 miRNAs and 12 mRNAs were randomly selected for qRT-PCR analysis. We compared the expression data of the four groups obtained by RNA-seq and qRT-PCR (Fig. 9, Fig. 10). The correlation between RNA-Seq (FPKM) and qPCR (2-ΔΔCt) results for the genes and miRNAs was calculated using log2-fold variation measurements. The expression change trend of most genes and miRNAs observed from qRT-PCR analysis was similar to that of high-throughput sequencing. Among the 12 genes, the correlation value (R2) of 10 genes was between 0.6362 and 0.991 (Fig. 9). For miRNAs, R2 of 10 miRNAs was ranging from 0.661 to 0.9894 (Fig. 10). To sum up, these data revealed that there was good consistency between the high-throughput sequencing and qRT-PCR, indicating the reliability and accuracy of small RNA and mRNA transcriptome sequencing.