Ginkgo (Ginkgo biloba L.), also known as Gongsunshu, is the oldest relict plant of extant gymnosperms dating from dinosaurs and is recognized as ‘living fossil’. Paleobotanists’ studies of fossil data have shown that seed plants experienced a decline in the early Permian, but the adaptable gymnosperms such as ginkgo and cedar developed slowly in the late Triassic and peaked in the Jurassic and early Cretaceous. Trichopitys, the hairy leaves during Autunian in southern France, are recognized as reliable ancestors of Ginkgo biloba, and the earliest and best preserved Ginkgo biloba is the Yima Ginkgo in Early Jurassic discovered in China (Tralau, 1968). In the early 1990s, morphological studies on male and female flower buds of Ginkgo biloba were reported. Then, the transcriptome sequencing and gene cloning of Ginkgo biloba flowering and secondary metabolites were reported in an endless stream (Zhang et al., 2013). Ginkgo biloba is a dioecious plant in sex differentiation. Its male reproductive organs are also called microsporangia, like the inflorescence of Rouyi, also known as male flowers (Zhang et al., 2001). The flower buds of leaves and male flowers on short branch are called ‘mixed male flower buds’, or called ‘male flower buds’ for short (Figure 1).

During plant growth, shoot apical meristem (SAM) continuously produces lateral branch organs. According to the morphological characteristics and functions of lateral branch organs, postembryonic development of plants can be divided into vegetative stage and reproductive stage. SAM produces leaves at vegetative stage and produces flowers at reproductive stage (Poethig, 2003). Flower transition, from vegetative stage to reproductive stage, is coordinated by multiple genetic pathways to respond to various developmental and environmental signals (Bäurle and Dean, 2006; Srikanth and Schmid, 2011). There are six genetic pathways involved in flowering regulation in A. thaliana, a model plant with relatively thorough research on flowering regulation. Among them, the research of photoperiod, temperature, vernalization, gibberellin and autonomous pathways is more thorough, while the research of age pathways remains to be further studied (Srikanth and Schmid, 2011). A. thaliana has a short growth cycle and is not strongly dependent on age pathways, but Ginkgo biloba has a very long life span and its childhood period is much longer than other plants, so it is particularly urgent to study the age-dependent pathway of flowering regulation. In this paper, transcriptome sequencing technology was used to conduct bioinformatics analysis of male flower buds in three age stages of Ginkgo biloba and screen differentially expressed genes, which can provide scientific basis for the study on age pathway of flowering regulation.

1 Results and Analysis

1.1 Sequencing data and quality control

A total of 57.45 Gb data were obtained by RNA-seq. The data amount of each sample reached 6.07 Gb, and the percentage of Q30 base was 94.59 % and above. The data of each sample was compared with the designated reference genome, and the comparison efficiency ranged from 85.64% to 90.19% (Table 1). The results showed that the sequencing quality is high and can be further analyzed.

1.2 Functional annotation of Ginkgo biloba Unigene

Unigenes sequences were compared with Nr (Nun-redundant Proteins) database, COG (Cluster of Orthologous Groups of Proteins) database, GO (Gene Ontology) database, KEGG (Kyoto Encyclopedia of Genes and Genomes) database, KOG (euKaryotic Orthologous Groups) database, Pfam protein database, Swiss-Prot (Swiss-Protein) database, and eggNOG (evolutionary genealogy of genes: Non-supervised Orthologous Groups) database to obtain annotation information (Table 2). Finally, in each database, a total of 35 058 genes obtained annotation information. Among the newly discovered 1 062 new genes, the number of genes with annotation information was 681. It can be seen that a large part of genes has not been annotated, which needs further exploration and research.

Ginkgo biloba Unigene was compared with the NR database to understand the evolutionary relationship between Ginkgo biloba and other species. According to the distribution of gene matching rate between Ginkgo biloba Unigene and other species in NR database (Table 3), it can be seen that Ginkgo biloba has the highest number of homologous genes with Picea sitchensis, which is also a gymnosperm. Apart from the high homology, it may also be related to the in-depth study of Picea sitchensis and the relatively clear annotations. In addition, 11 404 Unigenes were not compared with other species, which may be unique genes of Ginkgo biloba.

TIn organisms, different gene products construct a complex regulatory network, which can eventually perform biological functions through coordination with each other. Analyzing the annotations of gene pathway is helpful to understand the process of genes performing their functions. The annotation results of KEGG (Kyoto Encyclopedia of Genes and Genomes) showed (Table 4) that the Unigene of male flower buds of Ginkgo biloba involved 126 metabolic pathways in three age stages, among which starch and sucrose metabolism (304 genes), plant hormone signal transduction (302 genes), ubiquitin mediated proteolysis (222 genes) and plant circadian rhythm (72 genes) may play an indispensable role in flower bud induction. Other pathways are responsible for providing energy and raw materials and maintaining normal life activities during flower bud differentiation.

1.3 GO enrichment analysis of Ginkgo biloba Unigene

In order to further understand Ginkgo biloba Unigene enrichment locating on which biological function, pathway or cell, GO enrichment analysis was performed. The results showed that all Ginkgo biloba sequences were divided into three categories: Biological Process, Molecular Function, and Cellular Component, and that 35 990 entries, 27 079 entries and 15 534 entries were enriched, respectively. Among them, the Biological Process is rich in Biological Regulation, Developmental Process, Reproductive Process, Signaling, Growth, Reproduction and Rhythmic Process, which may be related to the age pathway of flowering regulation that this research focuses on. Take the first 50 items for drawing (Figure 2).

1.4 Screening of differentially expressed genes in age pathway of flowering regulation of Ginkgo biloba

The differential expression among sample groups was analyzed by edgeR (Figure 3). There are fewer up-regulated genes (37) and down-regulated genes (75) in junior male flower buds of Ginkgo biloba than in one-year male flower buds. While there are more up-regulated genes (592) and down-regulated genes (871) in junior male flower buds than in several-year male flower buds and more up-regulated genes (961) and down-regulated genes (1203) in one-year male flower buds than in several-year male flower buds.

KEGG pathway annotation and GO enrichment analysis were performed on differentially expressed genes, and 11 genes with significant expression differences that may be related to flowering regulation were further screened out (Table 5), and hierarchical clustering analysis was performed (Figure 4). After bioinformatics analysis, 11 differentially expressed genes were finally screened. The expression of SPLs (gene.Gb_23724, gene.Gb_03922) continued to decline during the transition of Ginkgo biloba from juvenile to adult, while AP2 (gene.Gb_00766) showed a continuous upward trend, which was consistent with the expression of miR156 and miR172. The expression patterns of the four genes with MADS-box domains (gene.Gb_01886, gene.Gb_15398, gene.Gb_28337, Gingko_newGene_2213) are not the same. SOC1 (gene.Gb_01886) is up-regulated and the rest is down-regulated. It indicated that multiple MADS-box genes may have redundant functions to regulate flowering in age pathway of flowering regulation of Ginkgo biloba. GA regulatory protein (gene.Gb_34467, gene.Gb_33214) is continuously down-regulated, and another GA regulatory protein (gene.Gb_28606) is up-regulated, indicating that in Ginkgo biloba, multiple regulatory proteins work together to regulate GA levels. The DELLA protein (gene.Gb_34644) regulated by GA signal showed a downward trend. It can be seen that GA-DELLAs played a certain regulatory role in the transition of Ginkgo biloba from non-flowering juvenile to flowering adult.

2 Discussion

In A. thaliana, SPLs with miR156 target sites are divided into two categories. One type encodes smaller proteins, such as SPL3, SPL4 and SPL5, which control flowering time and flowering period changes. The other type encodes larger proteins, including 8 genes such as SPL9 and SPL15, which promote leaf germination and flower formation (Cardon et al., 1997). In the stem end, both SPL3 and SPL9 can induce flowering by activating MADS-box genes including AP1, LFY, FUL and SOCl (Yamaguchi et al., 2009). And previous studies have shown that whether it is in short-day or long-day conditions, the overexpression of SPL3 or SPL9 of specific promoter in the stem end can lead to an early flowering phenotype (Wang et al., 2009).

AP2 transcriptional regulator is the other most important component of the age pathway. Studies have shown that AP2 can directly inhibit FT which up regulated by sunlight change (Castillejo and Pelaz, 2008). According to the difference of structural domains and binding sites, proteins containing AP2 domains can be divided into two types. One type is TEMPRANILLO1 (TEM1) and TEM2 proteins, which contain sites encoding AP2 domains and B3 type DNA binding domains, and have been shown to inhibit redundant flowering by inhibiting FT transcription (Castillejo and Pelaz, 2008). The other type is mainly targeting protein of miR172, including TARGET OF EAT (TOE) genes (TOE1 to TOE3), AP2, SCHLAFMU¨TZE (SMZ) and its homologous gene SCHNARCHZAPFEN (SNZ) (Aukerman and Sakai, 2003). AP2 can not only directly bind to SOC1 and AGAMOUS, the key genes for flowering regulation, and inhibit their transcription, but also can inhibit flowering and flower development by binding to miR172 gene locus and inhibiting its transcription (Yant et al., 2010). AP2 also acts as a dual-function transcription factor, directly binding to the genomic sites of the MADS-box family genes AGAMOUS-LIKE15 (AGL15) (Adamczyk et al., 2007) and miR156 (Wu et al., 2009), and directly activating other plant inhibitory factors, thereby regulating plant flowering.

Gibberellins (GA) is one of the most important plant hormones, which participates in many biological processes such as plant growth and development. It also plays an indispensable role in the flowering process of plants. Under short-day conditions, the addition of GA is a key factor to induce flowering of GA1 (GA REQUIRING 1) mutants (Wang, 2014). Under long-day conditions, GA also plays an important role in the flowering process (Wang et al., 2011). Between GA and miR156-SPL9, the florescence regulation is achieved by degrading the DELLAs protein (Dill et al., 2001). Under the specific promoters in leaf or stem end, the forced expression of GA-insensitive RGA or GA catabolism genes results in the late flowering phenotype (Galvão et al., 2012; Porri et al., 2012). The combination of RGA and SPL9 interferes with the activity of SPL9 transcriptome on the upstream of miR172, SOC1 and FUL. Therefore, DELLAs protein inhibits the phase transition process from growth to reproduction by inhibiting SOCl and FUL in the stem end. On the one hand, DELLA protein inhibits the expression of SOC1, transmits GA signals, and regulates flowering (Wilson et al., 1992). On the other hand, GA regulates MYB transcription factor by inhibiting the expression of miR159 through DELLA protein (Griffiths et al., 2006).

Because SOC1, FT, and LFY are regulated by multiple flowering pathways, they are also called flower integrons (Kardailsky et al., 1999; Kobayashi et al., 1999). As a downstream target gene of AP2, SOC1 is not only inhibited by AP2 in the age pathway, but also can be positively regulated by SPL9 (Samach et al., 2000). It can also receive signals from GA and is negatively regulated by the DELLA protein (Simpson and Dean, 2002).

The expression patterns of SPL3 (gene.Gb_23724) and AP2 (gene.Gb_00766) in the stem ends of Ginkgo biloba in three age stages are different from those of A. thaliana. The expression of GA receptor GID1 (GIBBERELLIC INSENSITIVE DWARF 1) is not significantly different in Ginkgo biloba, but GA regulatory protein shows significant up-regulation (gene.Gb_28606) and down-regulation (gene.Gb_34467, gene.Gb_33214) trends. The expression of DELLA protein RGA2 (gene.Gb_34644) is continuously down-regulated in Ginkgo biloba. The previous work of this study proved that GbMADS6 in Ginkgo biloba is a homologous gene of SOC1, which shows a continuous up-regulated expression pattern in three age stages of Ginkgo biloba, and promotes the flowering of Ginkgo biloba. In summary, it is indicated that the age-regulated flowering of Ginkgo biloba may be different from that of A. thaliana.

Based on transcriptome sequencing data, the key genes that are differentially expressed in three age stages of Ginkgo biloba was screened in this study. It can provide new data and theoretical basis for further cloning and functional verification of key genes, and for revealing the molecular mechanism of the flowering pathway in Ginkgo biloba and even other gymnosperms. It can also provide scientific basis for promoting early flowering, molecular breeding and spreading planting of Ginkgo biloba.

3 Materials and Methods

3.1 Materials

Juvenile male flower buds with 8 years’ growth cycle, one-year-old male flower buds with 21 years’ growth cycle and multi-year-old male flower buds with 45 years’ growth cycle in botanical garden in Central South University of Forestry and Technology were selected as experimental materials. The flower buds were collected at the beginning of the flower bud differentiation in early June. The lower xylem of the short branches was discarded, and the green part on the short branches, namely the male flower buds, was quickly frozen in liquid nitrogen and placed in the refrigerator at -80℃.

3.2 RNA extraction and library construction

The OMEGA kit was used to extract total RNA from male flower buds at three age stages of Ginkgo biloba, and the RNA concentration was measured by using NanoDrop 2000 (Thermo). The RNA Nano 6000 detection kit of the Agilent Bioanalyzer 2100 System (Agilent Technologies, California, USA) was used to evaluate RNA integrity. After the samples were qualified, NEBNext UltraTM RNA Library Prep Kit from Illumina (Nebraska, USA) was used to generate a sequencing library.

3.3 Transcriptome sequencing and quality control of sequencing data

Transcriptome sequencing was completed by Beijing Biomark Biotechnology Co., Ltd. After completing the sequencing, base quality score, base content distribution and data output statistics were used to control the probability of base calling errors, the presence or absence of AT and GC separation phenomenon and the percentage of Q30 bases in each sample.

3.4 Analysis of transcriptome data

The Clean Data obtained after splicing was performed sequence alignment with the specified reference genome, and then the obtained Mapped Data was performed quality control. Structural level analysis such as alternative splicing analysis, new gene mining and gene structure optimization was performed on Unigene. Based on the differential expression of genes, the analysis of differential expression of genes, function annotation of differential gene, functional enrichment and other expression level analysis were carried out.

3.5 Screening of differentially expressed genes

Choose edgeR Pvalue=0.01 and FC=2 as the screening criteria to draw the Venn diagram of differentially expressed genes in three age stages of Ginkgo biloba. On this basis, the functional annotation and enrichment analysis of the database of differentially expressed genes were carried out. The genes related to flowering were further selected, and the genes with the same or similar expression patterns were subjected to hierarchical cluster analysis.

Authors’ contributions

YTT was the experiment designer and executor of this study, and she wrote the first draft of the manuscript. FZ completed data analysis. YTT and WYQ are the creators and leaders of the project, directing data analysis, manuscript writing and revision. LM and other authors also gave a lot of help in this study. All authors read and approved the final manuscript.

Acknowledgements

This study was jointly funded by National Natural Science Foundation of China (31570682) and National Key R&D Program of China (2017YFD0600701).

References

Adamczyk B.J., Lehti-Shiu M.D., and Fernandez D.E., 2007, The MADS domain factors AGL15 and AGL18 act redundantly as repressors of the floral transition in Arabidopsis, The Plant J., 50(6): 1007-1019
https://doi.org/10.1111/j.1365-313X.2007.03105.x
PMid:17521410

Aukerman M.J., and Sakai H., 2003, Regulation of flowering time and floral organ identity by a microRNA and its APETALA2-like target genes, The Plant Cell, 15(11): 2730-2741
https://doi.org/10.1105/tpc.016238
PMid:14555699 PMCid:PMC280575

Bäurle I., and Dean C., 2006, The timing of developmental transitions in plants, Cell, 125(4): 655-664
https://doi.org/10.1016/j.cell.2006.05.005
PMid:16713560

Cardon G.H., Höhmann S., Nettesheim K., Saedler H., and Huijser P., 1997, Functional analysis of the Arabidopsis thaliana SBP‐box gene SPL3: a novel gene involved in the floral transition, The Plant J., 12(2): 367-377
https://doi.org/10.1046/j.1365-313X.1997.12020367.x
PMid:9301089

Castillejo C., and Pelaz S., 2008, The balance between CONSTANS and TEMPRANILLO activities determines FT expression to trigger flowering, Curr. Biol., 18(17): 367-377
https://doi.org/10.1046/j.1365-313X.1997.12020367.x
PMid:9301089

Dill A., Jung H.S., and Sun T., 2001, The DELLA motif is essential for gibberellin-induced degradation of RGA, P. Natl. A. Sci. India. B., 98(24): 14162-14167
https://doi.org/10.1073/pnas.251534098
PMid:11717468 PMCid:PMC61185

Galvão V.C., Horrer D., Küttner F., and Schmid M., 2012, Spatial control of flowering by DELLA proteins in Arabidopsis thaliana, Development, 139(21): 4072-4082
https://doi.org/10.1242/dev.080879
PMid:22992955

Griffiths J., Murase K., Rieu I., Zentella R., Zhang Z.L., Powers S.J., Gong F., Phillips A.L., Hedden P., Sun T.P., and Thomas S.G., 2006, Genetic characterization and functional analysis of the GID1 gibberellin receptors in Arabidopsis, The Plant Cell, 18(12): 3399-3414
https://doi.org/10.1105/tpc.106.047415
PMid:17194763 PMCid:PMC1785415

Kardailsky I., Shukla V.K., Ahn J.H., Dagenais N., Christensen S.K., Nguyen J.T., Chory J., Harrison M.J., Weigel D., 1999, Activation tagging of the floral inducer FT, Science, 286(5446): 1962-1965
https://doi.org/10.1126/science.286.5446.1962
PMid:10583961

Kobayashi Y., Kaya H., Goto K., Iwabuchi M., Araki T., 1999, A pair of related genes with antagonistic roles in mediating flowering signals, Science, 286(5446): 1960-1962
https://doi.org/10.1126/science.286.5446.1960
PMid:10583960

Poethig R.S., 2003, Phase change and the regulation of developmental timing in plants, Science, 301(5631): 334-336
https://doi.org/10.1126/science.1085328
PMid:12869752

Porri A., Torti S., Romera-Branchat M., and Coupland G., 2012, Spatially distinct regulatory roles for gibberellins in the promotion of flowering of Arabidopsis under long photoperiods, Development, 139(12): 2198-2209
https://doi.org/10.1242/dev.077164
PMid:22573618

Samach A., Onouchi H., Gold S.E., Ditta G.S., Schwarz-Sommer Z., Yanofsky M.F., and Coupland G., 2000, Distinct roles of CONSTANS target genes in reproductive development of Arabidopsis, Science, 288(5471): 1613-1616
https://doi.org/10.1126/science.288.5471.1613
PMid:10834834

Simpson G.G., and Dean C., 2002, Arabidopsis, the Rosetta stone of flowering time, Science, 296(5566): 285-289
https://doi.org/10.1126/science.296.5566.285
PMid:11951029

Srikanth A., and Schmid M., 2011, Regulation of flowering time: all roads lead to Rome, Cell Mol. Life Sci., 68(12): 2013-2037
https://doi.org/10.1007/s00018-011-0673-y
PMid:21611891

Tralau H., 1968, Evolutionary trends in the genus Ginkgo, Lethaia, 1: 63-101
https://doi.org/10.1111/j.1502-3931.1968.tb01728.x

Wang J.W., Czech B., and Weigel D., 2009, miR156-regulated SPL transcription factors define an endogenous flowering pathway in Arabidopsis thaliana, Cell, 138(4): 738-749
https://doi.org/10.1016/j.cell.2009.06.014
PMid:19703399

Wang J.W., 2014, Regulation of flowering time by the miR156-mediated age pathway, J. Exp. Bot., 65(17): 4723-4730
https://doi.org/10.1093/jxb/eru246
PMid:24958896

Wang R., Albani M.C., Vincent C., Bergonzi S., Luan M., Bai Y., Kiefer C., Castillo R., and Coupland G., 2011, Aa TFL1 confers an age-dependent response to vernalization in perennial Arabis alpina, The Plant Cell, 23(4): 1307-1321
https://doi.org/10.1105/tpc.111.083451
PMid:21498681 PMCid:PMC3101554

Wu G, Park M Y, Conway S R, Wang J W, Weigel D, and Poethig R S., 2009, The sequential action of miR156 and miR172 regulates developmental timing in Arabidopsis, Cell, 138(4): 750-759
https://doi.org/10.1016/j.cell.2009.06.031
PMid:19703400 PMCid:PMC2732587

Wilson R.N., Heckman J.W., and Somerville C.R., 1992, Gibberellin is required for flowering in Arabidopsis thaliana under short days, Plant Physiol., 100(1): 403-408
https://doi.org/10.1104/pp.100.1.403
PMid:16652976 PMCid:PMC1075565

Yamaguchi A., Wu M.F., Yang L., Wu G., Poethig R.S., and Wagner D., 2009, The microRNA-regulated SBP-Box transcription factor SPL3 is a direct upstream activator of LEAFY, FRUITFULL, and APETALA1, Dev. Cell., 17(2): 268-278
https://doi.org/10.1016/j.devcel.2009.06.007
PMid:19686687 PMCid:PMC2908246

Yant L., Mathieu J., Dinh T.T., Ott F., Lanz C., and Wollmann H., 2010, Orchestration of the floral transition and floral development in Arabidopsis by the bifunctional transcription factor APETALA2, The Plant Cell, 22(7): 2156-2170
https://doi.org/10.1105/tpc.110.075606
PMid:20675573 PMCid:PMC2929098

Zhang W.P., Shi J.K., Fan W.G., and Wen X.P., 2001, Morphological differentiation of Male flower bud of Ginkgo biloba L., Yuanyi Xuebao (Acta Horticulturae Sinica), 28(3): 255-258

Zhang H.J., Tan B.Y., and Cao F.L., 2013, Cloning and characterization of the geranylgeranyl pyrophosphate synthase gene from Ginkgo biloba Linn, Nanjing Linye Daxue Xuebao (Journal of Nanjing Forestry University ( Natural Sciences Edition)), 37(4): 8-12