Screening and identification of miRNAs related to sexual differentiation of strobili in Ginkgo biloba by integration analysis of small RNA, mRNA, and degradome sequencing

doi:10.21203/rs.3.rs-31410/v1

PDF

Research article

Screening and identification of miRNAs related to sexual differentiation of strobili in Ginkgo biloba by integration analysis of small RNA, mRNA, and degradome sequencing

https://doi.org/10.21203/rs.3.rs-31410/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 25 Aug, 2020

Read the published version in BMC Plant Biology →

You are reading this latest preprint version

Background: Ginkgo biloba, a typical dioecious plant, is a traditional medicinal plant widely planted. However, it has a long juvenile period, which severely affected the breeding and cultivation of superior ginkgo varieties.

Results: In order to clarify the complex mechanism of sexual differentiation in G. biloba strobili. Here, a total of 3,293 miRNAs were identified in buds and strobili of G. biloba, including 1,085 conserved miRNAs and 2,208 novel miRNAs using the three sequencing approaches of transcriptome, small RNA, and degradome. Comparative transcriptome analysis screened 4,346 and 7,087 differentially expressed genes (DEGs) in MB _vs_ FB and MS _vs_ OS, respectively. A total of 6,032 target genes were predicted for differentially expressed miRNA. The combined analysis of both small RNA and transcriptome datasets identiﬁed 51 miRNA-mRNA interaction pairs that may be involved in the process of G. biloba strobili sexual differentiation, of which 15 pairs were verified in the analysis of degradome sequencing.

Conclusions: The comprehensive analysis of the small RNA, RNA and degradome sequencing data in this study provided candidate genes and clarified the regulatory mechanism of sexual differentiation of G. biloba strobili from multiple perspectives.

Plant Physiology and Morphology

Plant Molecular Biology and Genetics

Differentially expressed genes

Ginkgo biloba

MiRNA

Strobili sexual determination

Target gene

Ginkgo biloba L., a perennial gymnosperm, is a living fossil in plant kingdom and called “Gonsun” tree in China. The extract of leaves and nuts from G. biloba contains a large number of active compounds, such as flavanol glycosides and terpene lactones, which can be used in the treatment of cardiovascular diseases [1]. At present, G. biloba is widely cultivated because it is commonly used in pharmacy, food, skin care products, landscape and other fields [2]. However, the composition of active ingredients in different varieties of G. biloba are very divers, which has a significant impact on the quality of G. biloba leaf extract. In addition, G. biloba is a typical dioecious plant with unique flowering features and has an important role in evolutionary history [3]. It has a long juvenile period and takes almost 30 years to blossom from seedling, which severely affected the breeding and cultivation of superior ginkgo varieties [4, 5]. Therefore, it is particularly important to study reproductive organ development and sexual differentiation of G. biloba.

The physiological and molecular mechanism of sexual differentiation and determination of plants have been extensively studied for many years. The mechanism of plant sexual determination mainly includes two systems: genetic sexual determination and environmental sexual determination. Dioecious plants are important materials for studying plant sexual determination and differentiation processes. Up to now, several genes involved in the process of plant sexual determination have been discovered and identified from dioecism species. Y-gene 1 (SLY1) was the first gene associated with sexual differentiation to be cloned and identified from Silene latifolia [6]. Subsequently, two MADS-box genes, SpAPETALA3 and SpPISTILLATA, were characterized to be involved in formation of male flower primordia in Spinacia oleracea [7]. In Cucumis sativus, CsACS2 gene was proved to inhibit male flower formation and promote female flower formation [8]. More recently, two important sex-differentiated genes, Shy Girl (SyGI) and Friendly boy (FrBy),, were cloned from Actinidia chinensis and identified to act for the maintenance of male functions [9, 10]. In addition, development of unisexual male flowers of Cucumis melo results from the expression of the sex determination gene, CmWIP1, which interacts with CmbZIP48 to inhibit carpel development [11].

Moreover, the MADS-box transcription factor is closely related to flower organ development and sexual differentiation. For example, five B-type MADS-box homologous genes involved in sexual differentiation were isolated from Silene alba but there is no relevant evidence to prove whether these five genes are located on the sexual chromosome [12]. In addition to key genes, most of phytohormones, such as auxin, ethylene (ETH), gibberellin (GA) and cytokinin (CTK), also played important role in sexual differentiation in plant [13]. ETH and CTK mainly promote female flower differentiation and have been widely used in fruit and vegetable industry [14]. IAA (indoleacetic acid) and auxin-like substances mainly act on the elongation of plant cells and promoting female phenotype. Spraying exogenous GA to C. sativus increased the male proportion in male and monoecious cucumber and induced the formation of male flowers in gynoecious plants [8]. Therefore, phytohormones play important roles in flower development and sexual differentiation in plant.

MicroRNA (miRNA) is a kind of endogenous non-coding single-stranded small RNA molecule composed of 18–24 nucleotides that is highly conservative in evolution [15]. Its main function is to regulate the expression level of target genes after transcription by cutting target genes or inhibiting the translation process of target genes [16]. In plants, miRNA participates in the regulation of many life processes, including the regulation of plant growth and development, morphogenesis, organ differentiation, tissue development, signal transduction, flowering and others, as well as the response to biotic and abiotic stresses [17, 18, 19]. For example, mutation of DCL1 reduces mature miRNA synthesis that changes the morphology of plant leaves, delays flowering, and leads to female sterility [20]. The growth and development of leaves are affected by miRNA396in Arabidopsis thaliana. Overexpressed miRNA396 resulted in changes in leaf size, decrease of stomatal density and enhancement of drought tolerance [21]. Furthermore, overexpression of miRNA159 led to downregulation of the target genes MYB33 and MYB65, resulting in delayed flowering and male sterility in plants [22, 23]. Also, miRNA plays an important role in plant sexual determination. MiRNA172 regulates the sexual differentiation of maize tassel by controlling the expression of the AP2 gene. In the hermaphroditic maize, miRNA172e targets the IDS1 transcription factor homologous to AP2, thereby preventing the development of pistils in male flowers and the development of intact male flowers. In female flowers, IDS1 transcription factors are not regulated by miRNA172e, so female flowers continue to form a normal pistil [24]. Recently, a Y-chromosome–encoded miRNA has been identified in Diospyros lotus, a dioecious plant,can regulate autosome-located WeGI gene, responsible for pollen abortion [25].

G. biloba is one of representative dioecious trees. Few studies provided information about genetic mechanism of flowering in G. biloba. Although our previous work cloned and function reported some genes involved in flower development including three MADS-box transcription factor genes [26, 27, 28] and two circadian-regulated genes (GbCO and GbFT) in photoperiodic pathway [4, 29]. However, little literature reported genetic regulation of sexual differentiation and sexual determination of G. biloba, especially for miRNA identification of strobili sexual differentiation and its regulatory mechanism. To further understand the mechanism of strobili sexual differentiation of G. biloba, we conducted small RNA sequencing, transcriptome sequencing, and degradation sequencing of male and female buds (MB and FB) as well as microstrobilus (MS) and ovulate strobilus (OS) in this study. Our data identified some important candidate miRNAs and their target genes that regulate sexual differentiationand provided novel insights into genetic mechanism of sexual differentiation of G. biloba.

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), GA₃ 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 GA₃ 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
Metabolic pathway /Development	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
ETH metabolic pathway	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 (R²) of 10 genes was between 0.6362 and 0.991 (Fig. 9). For miRNAs, R² 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.

Candidate miRNA-mRNA interactions participated in sex differentiation in G. biloba

Several literatures reported screening and characterization of a number of miRNA from G. biloba using high-throughput sequencing, which potential participated in seed germi seed germination [33], ovules development [3], primary metabolisms in leaves [34], and terpene lactone biosynthesis [35]. Nevertheless, little is known about the miRNAs involved in sex differentiation of G. biloba strobili to date. Recently, Du et al. [36] used transcriptome data analysis to screen some genes related to sex differentiation in Ginkgo. However, single transcriptome sequencing cannot fully illuminate the gene regulatory network involved in plant sex differentiation. The combined analysis of multiple sequencing could present itself as a versatile tool to analyze the mechanisms of plant life process. The complex mechanism of Ascochyta blight resistance in chickpea were analyzed by means of multi-group sequencing combination [37]. In Camellia, miR172-AP2 and miR156-SPLs regulatory pairs were identified as key regulatory nodes to promote the diversity of double flower forms by the multi-group analysis [38]. In this study, three high-throughput approaches, namely, transcriptomics, small RNA and degradome sequencings were used to elucidate the mechanisms of sexual differentiation in G. biloba. In addition, we identified miRNA-mRNA interaction pairs that may regulate the sex differentiation. Individual miRNA alone has no function and is easy to degrade. Only when miRNA combines with the Argonaute (AGO) protein to form the RNA-induced silencing complex (RISC) it can play an effective role, in which miRNA guides RISC to target-specific mRNA regions [39]. Plant miRNA can completely complement the coding regions of mRNAs, and inhibit the expression of genes by inducing the degradation of the complementing mRNAs [40]. A total of 51 miRNA-mRNA negative correlation regulatory relationships were identified in this research. The results of the present work indicated that miRNA might regulate sex differentiation by targeting transcription factor genes involved in the biosynthesis and signal transduction pathways of hormones, as well as flower development-related genes in G. biloba (Fig. 11).

Sexual differentiation of G. biloba might be affected by MiRNAs involved in ETH and GA biosynthesis pathways

ETH, one of plant hormones played important role in plant sexual determination, acts on the downstream of GA. ETH is a factor that promoted the development of female organs and inhibited the development of male organs [41]. Here, miR2950.1 and miR2950.2 targeted ethylene response transcription factor 2 (EMF2),, which was down-regulated in FB of G. biloba, as opposed to the results of other plants [41]. The differences in regulatory function in sex differentiation might be due to the different species studied.

MYB is a large family of transcription factors in plants and involved in many life activities such as plant secondary metabolism [42], hormone responses [43], leaf organ formation [44], flower development [45]. At present, some MYB transcription factors have been identified as regulators of flower development [46]. In A. thaliana, several MYBs were reported to be involved in jasmonic acid-mediated stamen maturation. Transcriptional analysis indicated that MYB108 and MYB24 have overlapping functions and acted downstream of MYB21 in a transcriptional cascade that mediated stamen and pollen maturation [47, 48]. Our findings was consistent with the results of previous studies, which demonstrated that GA could promote differentiation of male flowers [8]. The present study found that miR159.1, miR159.2, miR159.3, miR159.4 and miR858 were identified to target GAMYB1 in MB, FB, MS and OS. The GAMYB1 was down-regulated in FB and OS, indicating that the above-mentioned miRNAs target GAMYB1 gene to make GAMYB1 up-regulated in male strobili, thus promoting the maturation of G. biloba stamens. Moreover, it is well known that GA promotes the differentiation of the pistil. For example, maize, a monoecious plant, showed a phenotype in which the female flowers were transformed into male flowers or the flower primordium was biased toward pistil differentiation as a result of disruption of the GA signal transduction pathway or gene mutation during GA biosynthesis. These results also indicated endogenous GA promoted female flower differentiation in maize plants [41]. Likewise, our data also indicated that the PME13 gene participated in the biosynthesis of GA and was up-regulated in OS. In addition, the content of GA₃ in OS was significantly higher than that in MS, indicating that the expression level of this gene positively correlates with the GA₃ content in strobili of G. biloba. Therefore, these results suggested that GA was likely to play a role in producing female strobili in G. biloba.

Sexual differentiation might be affected by miRNAs functioned as regulator of IAA biosynthesis in G. biloba

IAA plays a key role in plant flower development such as flower primordia formation and flower organ differentiation [49]. In particular, most studies demonstrated that IAA promotes the expression of female flowers in plants [50]. For example, several ARFs and Aux/IAAs in IAA metabolic pathway have been identified to be essential for flower development in sugar apples [51]. In A. thaliana, ARF3 integrates the functions of AGAMOUS and APETALA2 in determining the flower meristem [52]. At the same time, two other members (ARF6 and ARF8) of ARF family were verified to play an important role in controlling the growth and development of plant vegetative organs and flower organs [53]. Furthermore, studies have shown that the loss of ARF18 activity causes the cell expansion in the silique wall, further promoting the development of siliques, increasing the photosynthetic products in siliques, and ultimately increasing the weight of the seeds [54]. In this study, six members of conserved miR160 family were identified to target ARF18.1 and ARF18.2 genes, which were differentially expressed only in MS _vs_ OS and were down-regulated in OS. And that, the expression level of miR160s is significantly negatively correlated with the expression level of the target gene identified in IAA biosynthesis among different organs (Table 2). Given that the IAA content in MS was higher than that in OS (Fig. 2), it can be speculated that miR160s participate in formation of male strobili through regulating IAA metabolism in G. biloba.

MiRNA regulating CTK signal transduction pathway affects sexual differentiation

Accumulated studies indicated that significant difference in CTK was usually detected between male and female flowers [55, 56, 57]. For example, treatment with exogenous CTKon Mercurialis annua L., a dioecious plant led the flower primordia to differentiate into pistil [55]. Similarly, treating exogenous CTK to the male flowers of spinach, hemp, grape and other plants caused changing of the male flowers into female flowers, indicating that CTK could promote the differentiation of flower primordia into female flowers [58]. One of the reason of effect CTK on plant sex differentiation could be regulatory function of sex-determining genes. For example, SyGI, the first sex-determining gene from kiwifruit, was found to inhibit gynoecial phenotype of through negatively regulating cytokinin signals [9]. In this study, miR2950.3 was identified to potentially regulate the signal transduction of CTK (Table 2). Moreover, miR2950.3 was differentially expressed only in MS _vs_ OS with up-regulation in OS and target AHP5 gene was reciprocally expressed. Given that the contents of ZT in the male strobili in initial and final stage was higher than that in the female strobili (Fig. 2), our data implied that miR2950.3 affected the metabolism of CTK, and then may promote the development of male strobili and inhibit the formation of female strobili in G. biloba.

MiRNAs involved in regulating flower development pathway affect sexual differentiation of G. biloba

Plants have no obvious secondary sexual characteristics, and the difference between male and female sex is mainly manifested in floral organs [59]. G. biloba is a typical dioecious plant, so genes related to flower development must participate in the process of sexual differentiation. MADS-box genes are transcription factors present in a high number and have conservative characteristics. They participate in almost all kinds of pathways to regulate plant development especially for flower development [60]. Most of the genes in the ABCDE model related to flower development also belong to MADS-box gene family. AP2, one of MADS-box transcription factor, was consider to function in regulation of the flowering time, regulation of the development of seeds and maintenance of the ability of plant apical meristem [61, 62]. In the case of MS _vs_ OS, two novel miRNAs (novel_miR_1944 and novel_miR_2329) simultaneously targeted the AP2 gene, which was up-regulated in OS. At the same time, the novel_miR_1858-ATR interaction pair was identified based on analysis of degradome sequencing. ATR wasup-regulated expressed in OS of G. biloba and was annotated to participate in flower development pathway. In addition, two other novel miRNAs (novel_miR_1881 and novel_miR_2912) were also found to be targeted genes related to flower development with up-regulated GIGANTEA and down-regulated POX53 in FB, respectively (Table 2). Taken together, our results inferred that the interaction between these miRNAs and target genes played a regulatory role in sexual differentiation of G. biloba. It can regulate sexual differentiation by regulating not only genes related to flower development directly, but also transcription factors indirectly.

This study is the first attempt to combine mRNA and miRNA expression data with degradome analysis to identify miRNAs involved in the regulation of sexual differentiation in G. biloba. Integration analysis of small RNA and mRNA were performed to identify 51 pairs of miRNA-mRNA interaction that may participate in regulating sexual differentiation in G. biloba. The finding was based on the screening of differentially expressed miRNA in combination with factors affecting sexual differentiation. The corresponding target mRNA also showed differential expression in male and female strobili of G. biloba. In total, 15 pairs among 51 miRNA-mRNA interaction pairs were confirmed by degradome sequencing. The comprehensive datasets of differential expression patterns of genes and miRNAs in male and female strobili presented a new insight to the sexual differential in G. biloba, and these sequencing data will also be a valuable resource for further study on G. biloba.

Plant material and RNA extraction

Thirty-two-year-old trees of G. biloba cv. ‘Jiafoshou’ is grown in the Ginkgo Garden of Yangtze University, China (30.35 °N, 112.14 °E). For the investigation of morphogenesis, the male and female strobili of G. biloba was observed and recorded with a digital camera from the March 8 to the April 12. The MB and FB harvested on March 13 as well as MS and OS harvested on April 5 were collected from G. biloba for extraction of total RNA to construct a sequencing library, and for qRT-PCR analysis and determinate the content of the hormone. A total nine independent Ginkgo trees were selected for the experiment, with 36 tree branches as the experimental block. Each biological replicate composed the samples selected from 12 independent branches from three ginkgo trees uniform in terms growth condition. A total of two biological repeats for the sequencing, and three biological repeats for the determination of the hormones content and the quantitative real-time PCR. The samples were collected, frozen in liquid nitrogen and stored in a refrigerator at –80 °C for further analysis. The total RNA was extracted to construct a sequencing library and perform qRT-PCR analysis using a PLANTpure kit (Aidlab, Beijing, China).

Scanning electron microscopy

The morphological of G. biloba cones were observed by SEM of the fresh MS and OS harvested on April 5. The ultrastructural analysis of MS and OS was performed according to the method described by Jin et al. [63] with following modification. The fresh sampleswere quickly fixed using 2.5% glutaraldehyde at 25 °C for 4 h and transferred to a refrigerator at 4°C for 12 h. The fixed samples were taken out, and rinsed 2–3 times with 0.1 M phosphate buffer and then were dehydrated in a graded ethanol series (30, 50, 70, 80, 90, and 100%, 15 min each step). The dehydrated materials were transferred to a clean beaker, and immersed under isoamyl acetate twice for 15 min each time. After that, the samples were placed into a critical point drying device (high-pressure closed container) with injecting liquid carbon dioxide at a critical state of 31 °C and at 72.8 standard atmosphere pressures. The samples were coated with a layer of gold and were observed under 15 kV accelerated voltage using QUANAT–200 scanning electron microscope.

Determination of phytohormone content

The extraction and content determination of IAA, GA₃ and ZT according to the method described by Onanuga et al. [64] with following modification. The preparation sample was as below: 5 g of fresh sample was ground into powder 4 °C, and then added 50 ml of precooled 80% methanol solution. The mixture of sample powder and methanol solution was stirred for 12h, ultrasonically crushed at 4 °C for 30 min, and then centrifuged at 12000 g/min for 10 min at 4 °C. Subsequently, the supernatant was extracted and concentrated to 10 mL using a rotary evaporator under reduced pressure at 35 °C; added equal volume petroleum ether for extraction for 3 times, and discard the organic phase. The aqueous extract were adjusted to pH 8.5 with 0.1 mol/L Na₂HPO₄ solution, and 0.2 g PVPP was added to shake well. After removing residue using a 0.45 μm filter membrane, the filtrate of aqueous extract was adjusted to pH 3.0 with 0.1 mol/L citric acid. Added equal volume of ethyl acetate for extraction for 3 times, discard the water phase, concentrate the remaining organic phase to about 2 mL using a rotary evaporator under reduced pressure at 35 °C, dissolved with methanol (chromatographic grade), and fixed the volume to 10 mL, then filtered with a 0.45 μm filter membrane to obtain the sample solution. High performance liquid chromatography (HPLC, Thermo Scientific™ UltiMate™ 3000, USA) with an Agilent C–18 chromatographic column (5 μm, 4.6 mm × 250 mm) was used for quantitative analysis of IAA, GA3 and ZT [64] with following modification. The mobile phase was methanol–1% acetic acid (40:60); flow rate: 0.6 mL/min; column temperature: 35 °C; injection volume: 2 μL; detection at 269 nm. Each sample was determined in six technical replicates.

Small RNA library preparation and sequencing

In total, eight RNA samples (MB–1, MB–2, FB–1, FB–2, MS–1, OS–1) were subjected to small RNA sequencing. The library was constructed using the Small RNA Library Prep Kit as following steps. T4 RNA Ligase 1 and T4 RNA Ligase 2 (truncated) were used to ligate the 3′-SR and 5′-SR Adaptor. Then, reverse transcription synthetic first chain and the PCR amplification were performed. PAGE gel was used to electrophoresis fragment screening purposes, rubber cutting recycling as the pieces get small RNA libraries. The library quality was assessed on the Agilent Bioanalyzer 2100 system. The raw reads obtained by sequencing contained adapter sequence or low-quality reads. In order to ensure the accuracy of information analysis, a series of quality control checks were carried out on the raw reads to obtain clean reads, including removing reads containing adapter, reads containing ploy-N and low-quality reads from raw data. The reads were trimmed and cleaned by removing the sequences smaller than 18 nt or longer than 30 nt. All the downstream analyses were based on clean data with high quality. The transcriptional level of small RNA was quantified according to Transcripts Per Kilobase of exon model per Million mapped reads (TPM) method [65]. The clean reads were compared with Silva database, GtRNAdb database, Rfam database and Repbase database, respectively, by using Bowtie software (v1.0.0, -v 0), and ncRNA such as ribosomal RNA (rRNA), transport RNA (tRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA) and repetitive sequences were filtered [66]. The remaining reads were used to detect conserved miRNA from miRbase (http://www.mirbase.org/), and predict novel miRNA by comparing with the reference Genomic seuqence of G. biloba [67]. The unique reads ranging from18–25 nt were collected and mapped to the reference genome by using the Bowtie software (v1.0.0, -v 0). All miRNA sequencing data are deposited in the NCBI SRA database (PRJNA590093).

Transcriptome library construction and screening of DEGs

Eight RNA samples same with samples subjected to small RNA sequencing was used for the mRNA transcriptome sequencing. The RNA-seq libraries were constructed using NEBNext Ultra^TM RNA Library Prep kit (NEB, E7530L, USA) according to the kit instruction. The raw reads were obtained by sequencing of Illumina HiSeq4000 platform. High quality sequences (clean reads) were obtained by eliminating low-quality reads and removing contaminated adapter and removing reads with an N ratio greater than 5% [68]. The clean reads were then mapped to the reference genome sequence by using Hisat2 tools soft to obtain the Mapped reads, which were further analyzed and annotated based on the reference genome sequence of G. biloba [67]. DEseq software (http://www.bioconductor.org/packages/release/bioc/html/DESeq.html) was applied to analyze differential expression gene (DEG) between sample groups [69]. Quantification of gene expression levels were estimated by fragments per kilobase of transcript per million fragments mapped (FPKM). The resulting P values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with an adjusted P-value < 0.05 found by DESeq were assigned as differentially expressed, and uses a Fold Change≥4 and FDR<0.001 as screening criteria to obtain differential expression gene sets between two samples. All transcriptome sequencing data are available in the NCBI SRA database (PRJNA590044).

Degradome library construction

Total RNA from the female and male strobili of G. biloba on March 13 and April 5 was mixed for the construction of the degradome library according to the kit instruction of NEBNext Ultra II RNA Library Prep Kit (NEB, E7775, USA). The degradome library was sequenced using an Illumina Hiseq 4000 platform [70]. Degradome reads were filtered using custom Perl script. Clean tags and cluster tags were obtained by removing adaptors and filtering low-quality raw tags. The cluster tags were compared to the reference genome to obtain the distribution of the tags on the genome [67]. Cluster tags and Rfam databases were compared to annotate non-coding RNA, and unannotated reads were used for degradation site analysis. The reads mapping to the sense strand of transcriptome were processed using the CleaveLand v4.4 pipeline to predict miRNA cleavage sites [32].

Identification of miRNAs

The conserved miRNAs were identified by comparing the reads on the reference genome with the conserved miRNA precursor in the miRBase database using the miRDeep2 software package (v2.0.5, -g –1 -l 250 -b 0) [71]. The possible precursors were obtained by comparing the position information of reads on the genome. Based on the distribution information on the genome of reads (miRNA production characteristics, mature, star, loop) and precursor structure energy information (randfold, v2.0, -s 99; RNAfold, v2.1.7, default), novel miRNAs were screened from the identiﬁed precursors using miRDeep-P, a probabilistic model-based miRNA prediction software especially designed for plant miRNAs [72].

Screening of differentially expressed miRNAs

A normalization factor calculated by DESeq R (v1.18.0, default) was used to normalize the raw read counts of miRNA followed by identiﬁcation of differentially expressed miRNA [73]. False Discovery Rate (FDR) and Fold Change (FC) are used as the criteria for screening differentially expressed miRNA, |log2 (FC)|> = 1; FDR< = 0.05.

Prediction of miRNA target genes and annotation of differentially expressed target genes

According to combination information of the conserved miRNAs and novel miRNAs, gene sequence of corresponding species, and high complementarity between miRNAs and their target genes, the prediction of target genes was performed using TargetFinder software (v1.0.0, -v 0) [74]. Using BLAST software (v2.2.26, -b 100 -v 100 -e 1e–5 -m 7 -a 2), the differentially expressed target genes were compared with the National Center for Biotechnology Information non-redundant (NCBI nr) [75], Swiss-Prot [76], Gene ontology (GO) [77], Kyoto Encyclopedia of Genes and Genomes (KEGG) [78] and Cluster of Orthologous Groups of proteins (COG) [79] databases for similarity (E-value 1e–05) to obtain annotation information.

qRT‑PCR validation of mRNA and miRNA

Forward specific primers were designed based on the sequence of mature miRNAs (Additional file 4: Table S4), and the reverse primers used the universal primers provided by the Mir-X miRNA First-Strand Synthesis Kit. The upstream and downstream primers (Tm: 60 °C) of the external gene U6 contained in the Mir-X miRNA First-Strand Synthesis Kit were used as the reference gene primers. The qRT-PCR reactions were performed using a TB GreenTM Advantage® kit (Takara, Dalian, China) in LineGene 9600 Plus (FQD–96A, China).

Specific primers of mRNA (sequences in Additional file 5: Table S2) were designed by software Primer 5.0 according to the sequences of the target genes. G. biloba 18S gene was used as the reference gene [80, 81]. The qRT-PCR reaction was executed according to the instruction manual of the AceQ® qPCR SYBR® Green Master Mix (Takara, Dalian). Experiments were performed with three independent biological triplicates and each triplicate consisted of three technical replicates. The expression of miRNAs and target gene mRNAs in flowers of G. biloba was calculated by the 2^-ΔΔCt method [82]. The standard error between biological replicates was analyzed by SPSS 22.0 (SPSS Inc., Chicago, IL, USA).

DEGs: differentially expressed genes; MB: male buds; FB: female buds; MS: microstrobilus; OS: ovulate strobilus; SLY1: Y-gene 1; SyGI: Shy Girl; FrBy: Friendly boy; IAA: indoleacetic acid; ETH: auxin, ethylene, GA: gibberellin; CTK: cytokinin; miRNA: MicroRNA; SEM: scanning electron microscope; ZT: Zeatin; ACS2: 1-aminocyclopropane–1-carboxylate synthase; CmWIP1: C2H2 zinc finger transcription factors named WIP proteins; bZIP48: basic region-leucine zipper; DCL1: DICER-LIKE1; GbCO: CONSTANS; GbFT: FLOWERING LOCUS T; qRT-PCR: quantitative real-time RCR; FPKM: fragments per kilobase of transcript per million fragments mapped; SPL: SQUMOSA PROMOTER BINDINGPROTEIN LIKE; AGO: Argonaute; RISC: RNA-induced silencing complex; EMF2: ethylene response transcription factor 2; AP2: APETALA2; PVPP: Polyvinylpolypyrrolidone; HPLC: High performance liquid chromatography; TPM: Transcripts Per Kilobase of exon model per Million mapped reads; rRNA: ribosomal RNA; tRNA: transport RNA; snRNA: small nuclear RNA; snoRNA: small nucleolar RNA; NCBI: National Center for Biotechnology Information; FDR: False Discovery Rate; FC: Fold Change; NCBI nr: National Center for Biotechnology Information non-redundant; GO: Gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; COG: Cluster of Orthologous Groups of proteins.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

The datasets generated and/or analysed during the current study are available in the NCBI repository, [NCBI SRA database (PRJNA590044) and NCBI SRA database (PRJNA590093)]. All data generated or analysed during this study are included in this published article Additional file. The datasets used and/or analyzed during the current study are available from the authors on reasonable request (Feng Xu, [email protected]; Xiao-Meng Liu, [email protected]).

Competing interests

The authors declare that they have no competing interests.

Funding

This study was supported by the National Natural Science Foundation of China (No. 31670608). The funding bodies did not play any role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Authors’ contributions

FX and SC conceived the idea and designed the experiment. XL, JY, ZC, and YL carried out various experiments. XL, WZ, and SK analyzed the data. XL and FX wrote the manuscript. All authors have read and approved the manuscript.

Acknowledgements

The authors would like to thank the head of the Ginkgo Garden of Yangtze University for us permission to collect all samples in this experiment.

Author information

¹College of Horticulture and Gardening, Yangtze University, Jingzhou 434025, Hubei, China. ²National R&D for Se-rich Agricultural Products Processing Technology, Wuhan Polytechnic University, Wuhan 430023, Hubei, China.³National Selenium Rich Product Quality Supervision and Inspection Center, Enshi 445000, Hubei, China. ⁴Research Institute for Special Plants, Chongqing University of Arts and Sciences 402160, ChongQing, China. ⁵Department of Agricultural Biotechnology, Seoul National University, Seoul 08826, Republic of Korea.

Ahlemeyer B, Krieglstein J. Neuroprotective effects of Ginkgo biloba extract. Cell Mol Life Sci. 2003;60(9):1779–1792.
Singh B, Kaur P, Sopichand GRD, Ahuja PS. Biology and chemistry of Ginkgo biloba. Fitoterapia. 2008;79(6):401–418.
Zhang Q, Li J, Sang Y, Xing S, Wu Q, Liu X. Identification and characterization of microRNAs in Ginkgo biloba var. epiphylla Mak. PLoS One. 2015a;10(5):e0127184.
Yan JP, Mao D, Liu XM, Wang LL, Xu F, Wang GY, et al. Isolation and functional characterization of a circadian-regulated CONSTANS homolog (GbCO) from Ginkgo biloba. Plant Cell Rep. 2017;36(9):1387–1399.
Mao D, Ye JB, Xu F. Advances of the Flowering Genes of Gymnosperms. Not Bot Horti Agrobo. 2019;47(1):1–9.
Delichere C, Veuskens J, Hernould M, Barbacar N, Mouras A, Negrutiu I, et al. SlY1, the first active gene cloned from a plant Y chromosome, encodes a WD-repeat protein. EMBO J. 1999;18(15):4169–4179.
Pfent C, Pobursky KJ, Sather DN, Golenberg EM. Characterization of SpAPETALA3 and SpPISTILLATA, B class floral identity genes in Spinacia oleracea, and their relationship to sexual dimorphism. Dev Genes Evol. 2005;215(3):132–142.
Zhang Y, Zhao GY, Li YS, Mo N, Zhang J, Liang Y. Transcriptomic analysis implies that GA regulates sex expression via ethylene-dependent and ethylene-independent pathways in cucumber (Cucumis sativus L.). Front Plant Sci. 2017;8:10.
Akagi T, Henry IM, Ohtani H, Morimoto T, Beppu K, Kataoka I, et al. A y-encoded suppressor of feminization arose via lineage-specific duplication of a cytokinin response regulator in kiwifruit. Plant Cell. 2018;30(4):780–795.
Akagi T, Pilkington SM, Varkonyi-Gasic E, Henry IM, Sugano SS, Sonoda M, et al. Two Y chromosome-encoded genes determine sex in kiwifruit. Nat Plants. 2019;5(8):801–809.
Eleblu JSY, Haraghi A, Mania B, Camps C, Rashid D, Morin H, et al. The gynoecious CmWIP1 transcription factor interacts with CmbZIP48 to inhibit carpel development. Sci Rep. 2019;9:15443.
Hardenack S, Ye D, Saedler H, Grant S. Comparison of MADS box gene expression in developing male and female flowers of the dioecious plant white campion. Plant Cell. 1994;6(12):1775–1787.
Khryanin VN. Role of phytohormones in sex differentiation in plants. Russ J Plant Physl. 2002;49(4):545–551.
Gerashchenkov GA, Rozhnova NA. The involvement of phytohormones in the plant sex regulation. Russ J Plant Physl. 2013;60(5):597–610.
Park MY, Wu G, Gonzalez-Sulser A, Vaucheret H, Poethig RS. Nuclear processing and export of microRNAs in Arabidopsis. Proc Natl Acad Sci USA. 2005;102(10):3691–3696.
Barrel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116(2):281–297.
Guo HS, Xie Q, Fei JF, Chua NH. MicroRNA directs mRNA cleavage of the transcription factor NAC1 to downregulate auxin signals for Arabidopsis lateral root development. Plant Cell. 2005;17(5):1376–1386.
Medina C, da Rocha M, Magliano M, Ratpopoulo A, Revel B, Marteu N, et al. Characterization of microRNAs from Arabidopsis galls highlights a role for miR159 in the plant response to the root-knot nematode Meloidogyne incognita. New Phytol. 2017;216(3):882–896.
Liu S, Mi X, Zhang R, An Y, Zhou Q, Yang T, et al. Integrated analysis of miRNAs and their targets reveals that miR319c/TCP2 regulates apical bud burst in tea plant (Camellia sinensis). Planta. 2019;250(4):1111–1129.
Schauer SE, Jacobsen SE, Meinke DW, Ray A. DICER-LIKE1: blind men and elephants in Arabidopsis development. Trends Plant Sci. 2002;7(11):487–491.
Liu DM, Song Y, Chen ZX, Yu DQ. Ectopic expression of miR396 suppresses GRF target gene expression and alters leaf growth in Arabidopsis. Physiol Plantarum. 2009;136(2):223–236.
Millar AA, Waterhouse PM. Plant and animal microRNAs: similarities and differences. Funct Integr Genomic. 2005;5(3):129–135.
Palatnik JF, Allen E, Wu XL, Schommer C, Schwab R, Carrington JC, Weigel D. Control of leaf morphogenesis by microRNAs. Nature. 2003;425(6955):257–263.
Banks JA. MicroRNA, sex determination and floral meristem determinacy in maize. Genome Biol. 2008;9(1):204.
Akagi T, Henry IM, Tao R, Comai L. SPL. Science. 2014;346(6209):646–650.
Cheng S, Cheng J, Xu F, Ye J, Wang XH. Molecular cloning and expression analysis of a putative e class MADS-box gene, GbSEP, from Ginkgo biloba. J Anim Plant Sci. 2016;26(1):253–260.
Wang XH, Cheng JH, Xu F, Li XX, Zhang WW, Liao YL, et al. Molecular cloning and expression analysis of a MADS-Box gene (GbMADS2) from Ginkgo biloba. Not Bot Horti Agrobo. 2015b:43(1):19–24.
Yang F, Xu F, Wang XH, Liao YL, Chen QW, Meng XX. Characterization and functional analysis of a MADS-box transcription factor gene (GbMADS9) from Ginkgo biloba. Sci Hortic. 2016;212:104–114.
Wang LL, Yan JP, Zhou X, Cheng SY, Chen ZX, Song QL, et al. GbFT, a FLOWERING LOCUS T homolog from Ginkgo biloba, promotes flowering in transgenic Arabidopsis. Sci Hortic. 2019;247:205–215.
Wang LW, Liu HH, Li DT, Chen HB. Identification and characterization of maize microRNAs involved in the very early stage of seed germination. BMC Genomics. 2011;12:154.
Aryal R, Ming R. Sex determination in flowering plants: Papaya as a model system. Plant Sci. 2014;217:56–62.
Addo-Quaye C, Miller W, Axtell MJ. CleaveLand: a pipeline for using degradome data to find cleaved small RNA targets. Bioinformatics. 2009;25(1):130–131.
Cui J, Zhao J, Zhao J, Xu H, Wang L, Jin B. Cytological and miRNA expression changes during the vascular cambial transition from the dormant stage to the active stage in Ginkgo biloba L. Trees. 2016;30(6):2177–2188.
Wang L, Zhao JG, Zhang M, Li WX, Luo KG, Lu ZG, et al. Identification and characterization of microRNA expression in Ginkgo biloba L. leaves. Tree Genet Genoms. 2015a;11(4):76.
Ye J, Zhang X, Tan J, Xu F, Cheng S, Chen Z, et al. Global identification of Ginkgo biloba microRNAs and insight into their role in metabolism regulatory network of terpene trilactones by high-throughput sequencing and degradome analysis. Ind Crop Prod. 2020;148:112289.
Du SH, Sang YL, Liu XJ, Xing SY, Li JH, Tang HX, et al. Transcriptome profile analysis from different sex types of Ginkgo biloba L. Front Plant Sci. 2016;7:871.
Garg V, Khan AW, Kudapa H, Kale SM, Chitikineni A, Sun QW, et al. Integrated transcriptome, small RNA and degradome sequencing approaches provide insights into Ascochyta blight resistance in chickpea. Plant Biotechnol J. 2019;17(5):914–931.
Li X, Li J, Fan Z, Liu Z, Tanaka T, Yin H. Global gene expression defines faded whorl specification of double flower domestication in Camellia. Sci Rep. 2017;7:3197.
Schwarz DS, Hutvagner G, Du T, Xu ZS, Aronin N, Zamore PD. Asymmetry in the assembly of the RNAi enzyme complex. Cell. 2003;115(2):199–208.
Yu B, Wang H. Translational inhibition by microRNAs in plants. Prog Mol Subcell Biol. 2010;50:41–57.
Tanurdzic M, Banks JA. Sex-determining mechanisms in land plants. Plant Cell. 2004;16:S61-S71.
Lea US, Slimestad R, Smedvig P, Lillo C. Nitrogen deficiency enhances expression of specific MYB and bHLH transcription factors and accumulation of end products in the flavonoid pathway. Planta. 2007;225(5):1245–1253.
Zhang TQ, Zhao YL, Wang YC, Liu ZY, Gao CQ. Comprehensive analysis of MYB gene family and their expressions under abiotic stresses and hormone treatments in Tamarix hispida. Front Plant Sci. 2018;9:1303.
Gu ZY, Zhu J, Hao Q, Yuan YW, Duan YW, Men SQ, et al. A novel R2R3-MYB transcription factor contributes to petal blotch formation by regulating organ-specific expression of pschs in tree peony (Paeonia suffruticosa). Plant Cell Physiol. 2019;60(3):599–611.
Saha G, Park JI, Ahmed NU, Kayum MA, Kang KK, Nou IS. Characterization and expression profiling of MYB transcription factors against stresses and during male organ development in Chinese cabbage (Brassica rapa ssp pekinensis). Plant Physiol Bioch. 2016;104:200–215.
Vimolmangkang S, Han YP, Wei GC, Korban SS. An apple MYB transcription factor, MdMYB3, is involved in regulation of anthocyanin biosynthesis and flower development. BMC Plant Biol. 2013;13:176.
Mandaokar A, Browse MJ. MYB108 acts together with MYB 24 to regulate jasmonate-mediated stamen maturation in Arabidopsis. Plant Physiol. 2009;149(2):851–862.
Song SS, Qi TC, Huang H, Ren QC, Wu DW, Chang CQ, et al. The Jasmonate-ZIM domain proteins interact with the R2R3-MYB transcription factors MYB21 and MYB24 to affect Jasmonate-Regulated stamen development in Arabidopsis. Plant Cell. 2011;23(3):1000–1013.
Alabadi D, Blazquez MA, Carbonell J, Ferrandiz C, Perez-Amador MA. Instructive roles for hormones in plant development. Int J Dev Bio. 2009;53(8–10):1597–1608.
Xu G, Huang J, Yang Y, Yao YA. Transcriptome analysis of flower sex differentiation in Jatropha curcas L. using RNA sequencing. PLoS One. 2016;11(2):e0145613.
Liu KD, Feng SX, Pan YL, Zhong JD, Chen Y, Yuan CC. Transcriptome analysis and identification of genes associated with floral transition and flower development in sugar apple (Annona squamosa L.). Front Plant Sci. 2016;7:1695.
Liu XG, Dinh TT, Li DM, Shi BH, Li YP, Cao XW, et al. AUXIN RESPONSE FACTOR3 integrates the functions of AGAMOUS and APETALA2 in floral meristem determinacy. Plant J. 2014b;80(4):629–641.
Liu N, Wu S, Van Houten J, Wang Y, Ding B, Fei ZJ, et al. Down-regulation of AUXIN RESPONSE FACTORS 6 and 8 by microRNA 167 leads to floral development defects and female sterility in tomato. J Exp Bot. 2014a;65(9):2507–2520.
Liu J, Hua W, Hu Z, Yang H, Zhang L, Li R, et al. Natural variation in ARF18 gene simultaneously affects seed weight and silique length in polyploid rapeseed. Proc Natl Acad Sci USA. 2015;112(37):E5123-E5132.
Dellaporta SL, Calderon-Urrea DA. Sex determination in flowering plants. Plant Cell. 1993;5(10):1241–1251.
Hui WK, Yang YT, Wu GJ, Peng CC, Chen XY, Zayed MZ. Transcriptome profile analysis reveals the regulation mechanism of floral sex differentiation in Jatropha curcas L. Sci Rep. 2017;7:16421.
Ni J, Shah FA, Liu WB, Wang QJ, Wang DD, Zhao WW, et al. Comparative transcriptome analysis reveals the regulatory networks of cytokinin in promoting the floral feminization in the oil plant Sapium sebiferum. BMC Plant Biol. 2018;18:96.
Cheng ZJ, Zhu SS, Gao XQ, Zhang XS. Cytokinin and auxin regulates WUS induction and inflorescence regeneration in vitro in Arabidopsis. Plant Cell Rep. 2010;29(8):927–933.
Li DD, Sheng YY, Niu HH, Li Z. Gene interactions regulating sex determination in cucurbits. Front Plant Sci. 2019;10:1231.
Fujiwara S, Nakagawa M, Oda A, Kato K, Mizoguchi T. Photoperiod pathway regulates expression of MAF5 and FLC that encode MADS-box transcription factors of the FLC family in Arabidopsis. Plant Biotechnol. 2010;27(5):447–454.
Wurschum T, Gross-Hardt R, Laux T. APETALA2 regulates the stem cell niche in the Arabidopsis shoot meristem. Plant Cell. 2006;18(2):295–307.
Yant L, Mathieu J, Dinh TT, Ott F, Lanz C, Wollmann H, et al. Orchestration of the floral transition and floral development in Arabidopsis by the bifunctional transcription factor APETALA2. Plant Cell. 2010;22(7):2156–2170.
Jin B, Wang L, Wang J, Teng NJ, He XD, Mu XJ, Wang YL. The structure and roles of sterile flowers in Viburnum macrocephalum f. keteleeri (Adoxaceae). Plant Biol. 2010;12(6):853–862.
Onanuga AO, Jiang P, Adl S. Determination of Endogenous Hormones Content in Cotton Varieties (Gossypiumhirsutum) as Influenced by Phosphorus and Potassium Nutrition. J Agr Sci. 2012;4:7.
Fahlgren N, Howell MD, Kasschau KD, Chapman EJ, Sullivan CM, Cumbie JS, et al. High-throughput sequencing of Arabidopsis microRNAs: evidence for frequent birth and death of MIRNA genes. PLoS One. 2007;2(2):e219.
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(3):R25.
Guan R, Zhao Y, Zhang H, Fan GY, Liu X, Zhou WB, et al. Draft genome of the living fossil Ginkgo biloba. GigaScience. 2016;5:49.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a ﬂexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–2120.
Wang LK, Feng ZX, Wang X, Wang XW, Zhang XG. Degseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26(1):136–138.
German MA, Pillay M, Jeong DH, Hetawal A, Luo SJ, Janardhanan P, et al. Global identiﬁcation of microRNA-target RNA pairs by parallel analysis of RNA ends. Nat Biotechnol. 2008;26(8):941–946.
Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. MirDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40(1):37–52.
Yang XZ, Li L. MiRDeep-P: a computational tool for analyzing the microRNA transcriptome in plants. Bioinformatics. 2011;27(18):2614–2615.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Allen E, Xie Z, Gustafson AM, Carrington JC. MicroRNA-directed phasing during trans-acting siRNA biogenesis in plants. Cell. 2005;121(2):207–221.
Deng YY, Li JQ, Wu SF, Zhu YP, Chen YW, He FC. Integrated NR database in protein annotation system and its localization. Comput En. 2006;32:71–74.
Apweiler R, Bairoch A, Wu CH, Barker WC, Boeckmann B, Ferro S, et al. UniProt: the universal protein knowledgebase. Nucleic Acids Res. 2004;32:D115-D119.
Ashburner M, Ball CA, Blake JA, Botstein D, Cherry JM. Gene Ontology: tool for the unification of biology. Nat Genet. 2000;25(1):25–29.
Kanehisa M, Goto S, Kawashima S, Okuno Y, Hattori M. The KEGG resource for deciphering the genome. Nucleic Acids Res. 2004;32:D277-D280.
Tatusov RL, Galperin MY, Natale DA, Koonin EV. The COG database: a tool for genome-scale analysis of protein functions and evolution. Nucleic Acids Res. 2000;28(1):33–36.
Cheng S, Yan J, Meng X, Zhang W, Liao Y, Ye J, et al. Characterization and expression patterns of a cinnamate–4-hydroxylase gene involved in lignin biosynthesis and in response to various stresses and hormonal treatments in Ginkgo biloba. Acta Physiol Plant. 2018;40(1):1–15.
Ye J, Cheng S, Zhou X, Chen ZX, Kim SU, Tan JP, et al. A global survey of full-length transcriptome of Ginkgo biloba reveals transcript variants involved in flavonoid biosynthesis. Ind Crop Prod. 2019;139:111547.
Schmittgen TD, Livak KJ. Analyzing real-time PCR data by the comparative CT method. Nat Protoc. 2008;3(6):1101–1108.

PDF

Journal Publication

published 25 Aug, 2020

Read the published version in BMC Plant Biology →

Submission checks completed at journal
10 Jun, 2020
Editor invited by journal
10 Jun, 2020

You are reading this latest preprint version

Screening and identification of miRNAs related to sexual differentiation of strobili in Ginkgo biloba by integration analysis of small RNA, mRNA, and degradome sequencing

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Results

Morphological changes of strobili during the flower development

Comparison of phytohormone contents of female and male strobili of G. biloba

Overview of small RNA deep-sequencing data

Differential expression of miRNAs

Screening of differentially expressed genes (DEGs)

Regulatory analysis between miRNAs and mRNA

Identification of differentially expressed miRNAs and target genes (DEGs) involved in regulating sexual differentiation in G. biloba

Discussion

Candidate miRNA-mRNA interactions participated in sex differentiation in G. biloba

Sexual differentiation of G. biloba might be affected by MiRNAs involved in ETH and GA biosynthesis pathways

MiRNA regulating CTK signal transduction pathway affects sexual differentiation

MiRNAs involved in regulating flower development pathway affect sexual differentiation of G. biloba

Conclusion

Methods

Plant material and RNA extraction

Scanning electron microscopy

Determination of phytohormone content

Identification of miRNAs

Screening of differentially expressed miRNAs

Prediction of miRNA target genes and annotation of differentially expressed target genes

qRT‑PCR validation of mRNA and miRNA

Abbreviations

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1