Temporal Dynamic Transcriptome Landscape Reveals Regulatory Network During the Early Differentiation of Female Strobilus Buds in Ginkgo biloba

Reproductive bud differentiation is one of the most critical events for the reproductive success of seed plants. Yet, our understanding of genetic basis remains limited for the development of the reproductive organ of gymnosperms, namely, unisexual strobilus or cone, leaving its regulatory network largely unknown for strobilus bud differentiation. In this study, we analyzed the temporal dynamic landscapes of genes, long non-coding RNAs (lncRNAs), and microRNAs (miRNAs) during the early differentiation of female strobilus buds in Ginkgo biloba based on the whole transcriptome sequencing. Results suggested that the functions of three genes, i.e., Gb_19790 (GbFT), Gb_13989 (GinNdly), and Gb_16301 (AG), were conserved in both angiosperms and gymnosperms at the initial differentiation stage. The expression of genes, lncRNAs, and miRNAs underwent substantial changes from the initial differentiation to the enlargement of ovule stalk primordia. Besides protein-coding genes, 364 lncRNAs and 15 miRNAs were determined to be functional. Moreover, a competing endogenous RNA (ceRNA) network comprising 10,248 lncRNA-miRNA-mRNA pairs was identified, which was highly correlated with the development of ovulate stalk primordia. Using the living fossil ginkgo as the study system, this study not only reveals the expression patterns of genes related to flowering but also provides novel insights into the regulatory networks of lncRNAs and miRNAs, especially the ceRNA network, paving the way for future studies concerning the underlying regulation mechanisms of strobilus bud differentiation.

Introduction

A flower represents a unique feature of angiosperms, which is absent in its sister group, gymnosperms. Instead, the reproductive organ of gymnosperms, namely, the strobilus or cone, lacks the basic structure found in angiosperms and is more primitive than that of angiosperms. Compared with angiosperms, which are mostly bisexual, unisexual strobilus (dioecious or monoecious) is common among extant gymnosperm species. Of note, approximately 65% of gymnosperms are dioecious and 1% are monoecious (Walas et al., 2018), while only 6% of angiosperms are dioecious and 7% are monoecious (Diggle et al., 2011). The research focused on the differentiation of strobilus buds will lead to a comprehensive understanding of floral evolution as well as the determination and differentiation of unisexual reproductive organs.

Reproductive bud differentiation is one of the key stages in the life cycle of seed plants, which is regulated by both the intrinsic and the extrinsic factors (Fornara et al., 2010; Blumel et al., 2015; Cho et al., 2017). It comprises three phases, namely, the vegetative-to-reproductive transition phase, the organ patterning phase, and the organ development phase (Luo et al., 2013). The regulatory network of flower bud differentiation has been extensively studied in some model systems of angiosperms. The floral transition from the vegetative growth, i.e., floral induction, is mainly involved with six pathways in Arabidopsis thaliana, i.e., the vernalization pathway, the photoperiod pathway, the ambient temperature pathway, the age pathway, the autonomous pathway, and the gibberellin pathway (Fornara et al., 2010). The downstream integrators of the six upstream pathways rapidly promote floral development and activate meristem identity genes (Weigel et al., 1992; Moon et al., 2003; Navarro et al., 2011). Generally, A+E, A+B+E, B+C+E, C+E, and C+D+E control the formation of sepal, petal, stamen, carpel, and ovule, respectively (Theissen and Saedler, 2001). Besides protein-coding genes, two key microRNAs (miRNAs) (miR156 and miR172) acted as count-down and count-up timers in the floral induction, respectively (Aukerman and Sakai, 2003; Wu and Poethig, 2006; Chuck et al., 2007; Zhang et al., 2015). COOLAIR is the first well-characterized long non-coding RNA (lncRNA) that regulates vernalization by silencing FLOWERING LOCUS C (FLC) transcription transiently in Arabidopsis (Swiezewski et al., 2009). Recently, the competing endogenous RNA (ceRNA) hypothesis has gained increasing attention, which shows lncRNAs, pseudogenes, and circular RNAs (circRNAs) can act as ceRNAs to compete for binding miRNAs through miRNA response elements (MERs) and decrease the suppression of miRNAs to their target mRNA (Salmena et al., 2011). Two lncRNAs, i.e., XLOC_0063639 and XLOC_007072, were found to sequester miR160 or miR164 as ceRNAs in the sexual reproduction of rice (Zhang et al., 2014). A ceRNA network containing 78 lncRNAs, 397 miRNAs, and 32 mRNA was involved in the regulation of tomato flowering (Yang et al., 2019).

Compared with angiosperms, the genetic basis of strobilus bud differentiation in gymnosperms has been paid less attention. Nevertheless, previous studies have reported a series of homologous genes related to the differentiation of strobilus buds in gymnosperms. Homologous genes of LEAFY (LFY) have double-copy in gymnosperms, which suggested that they have undergone functional differentiation (Frohlich and Parker, 2000; Shindo et al., 2001; Dornelas and Rodriguez, 2005; Guo et al., 2005). SUPPRESSOR OF OVEREXPRESSION OF CO 1-like (SOC1-like) genes have been cloned in some gymnosperms, such as Picea abies, Gnetum gnemon, and Pinus radiata (Tandre et al., 1995; Walden et al., 1998; Winter et al., 1999). The BC/D model has been established in gymnosperms. The homeotic C and D class genes are used to distinguish reproductive organs from vegetative organs, while B class genes play a key role in specifying the identity of microsporophyll (Theissen et al., 2000). Yet, the regulatory network of strobilus bud differentiation remains unknown (Mao et al., 2019).

The living fossil ginkgo (Ginkgo biloba L.) is the sole living species of Ginkgophytes in gymnosperms whose female reproductive organ, i.e., the ovulate strobilus or the female cone, represents one of the most primitive compound strobilus types, making it of great importance in understanding the evolution of strobili and flowers (Fu and Yang, 1993). The differentiation process of gingko female strobilus buds (FSBs) consists of six stages, i.e., the vegetative growth, the initial differentiation, the exuberant differentiation, the integument differentiation, the nucellus differentiation, and the collar differentiation stages (Shi et al., 1998; Douglas et al., 2007). A total of 29 genes were determined to be closely related to three differentiation stages during the differentiation process (He et al., 2018). The homologous genes of FLOWERING LOCUS T (FT), SOC1, and LFY were cloned, respectively (Zhang et al., 2002; Guo et al., 2005; Wang et al., 2019; Feng et al., 2021). Of note, 11 and 37 MIKC^c-type MADS-box genes were identified from a draft ginkgo genome and a new one, respectively (Guan et al., 2016; Chen et al., 2017; Liu et al., 2021). However, the functions of the homologous genes mentioned above have not been fully verified, and no attempts have been made to investigate the regulatory network of those genes or to pinpoint the role of non-coding RNAs (ncRNAs) during the differentiation of gingko FSBs.

In this study, we identified genes, lncRNAs, and miRNAs associated with the early differentiation of ginkgo FSBs using a whole-transcriptome sequencing strategy. We aimed (1) to explore the expression patterns of homologous genes (such as FT, SOC1, and LFY) related to the differentiation of ginkgo FSBs, (2) to identify the functional lncRNAs and miRNAs and draw the temporal dynamic expression landscapes, and (3) to explore the regulatory functions of ceRNA networks if exist. Our results provide new insights for further investigations on the regulatory networks of strobilus buds differentiation in gymnosperms, a land plant group with great economical, ecological, and evolutionary importance.

Materials and Methods

Plant Materials and Morphological Examination

Ginkgo FSBs were sampled from May to August in 2020 and 2021 from female trees with similar ages and growth conditions at Zhejiang University (Hangzhou, China) based on long-term observation. The fresh samples were fixed using 2.5% glutaraldehyde with the phosphate buffer (0.1 M, pH = 7.0). Post-fixation with 1% OsO₄ in the phosphate buffer for 1–2 h precedes three times of washing for 15 min. Samples were subjected to gradient dehydration using ethanol solutions (30, 50, 70, 80, 90, and 95%, respectively) for 15 min and 100% ethanol twice for 20 min. Then, samples were dehydrated in a Hitachi Model HCP-2 critical point dryer (Hitachi, Tokyo, Japan) before being embedded with gold-palladium for 4–5 min in Hitachi E-1010 ion sputter (Hitachi, Tokyo, Japan). Later, we observed those samples using a Hitachi SU-8010 Scanning Electron Microscope (SEM) (Hitachi, Tokyo, Japan).

Based on the SEM observations, FSB samples for RNA sequencing were collected at three stages of strobilus buds each from three female trees according to the number and characteristics of ovulate stalk primordia in the FSBs (Shi et al., 1998; Douglas et al., 2007), i.e., the vegetative growth (ST1, sampled on May 6), the initial differentiation (ST2, sampled on May 26), and the exuberant differentiation (ST3). Given that the exuberant differentiation stage lasted for a long period, two samples with large morphological differences were sampled in each tree (ST3-1, sampled on June 19; ST3-2, sampled on July 28).

Total RNA Extraction, Library Preparation, RNA Sequencing, and Quality Control

Total RNA was extracted using TRIzol Reagent (Invitrogen, Waltham, MA, United States) according to the manufacturer’s instructions. Subsequently, total RNA was qualified and quantified using NanoDrop ND-1000 (NanoDrop, Wilmington, DE, United States) and Agilent 2011 bioanalyzer (Agilent Technologies, Santa Clara, CA, United States).

For strand-specific libraries, approximately 1 μg of total RNA per sample was treated using the Ribo-Zero Plant Kit (Illumina, San Diego, CA, United States) to remove rRNA. The strand-specific libraries were constructed following RNA fragmentation, 1st cDNA strand synthesis, 2nd cDNA strand synthesis, end repair, tailing, and adaptor ligation. Finally, the qualified libraries were sequenced paired end on the Hiseq 4000 platform at BGI (Shenzhen, China). For small RNA libraries, approximately 1 μg of total RNA was used for each sample, and small RNA sequencing (sRNA-seq) was performed on the BGISEQ-500 platform (BGI, Shenzhen, China).

We removed adaptors, low quality, and ambiguous bases of raw data from strand-specific RNA-seq (ssRNA-seq) libraries and sRNA-seq libraries using SOAPnuke (Chen et al., 2018). Then, we discarded sRNA-seq reads shorter than 18 nt. After quality control using FastQC¹, clean reads were used for the subsequent analyses.

Reads Mapping and Transcript Quantification

The clean reads from ssRNA-seq were mapped to the updated draft ginkgo genome (Guan et al., 2019) using HISAT2 (Kim et al., 2019) with default parameters except “–rna-strandness RF.” StringTie (Pertea et al., 2015) was used to estimate the expression levels of all transcripts based on annotations of the ginkgo genome. The expression level of each gene was normalized to fragments per kilobase of transcript per million mapped reads (FPKM).

Identification of lncRNAs

The mapped reads from ssRNA-seq were assembled using StringTie. Then, all StringTie assemblies were merged into a single transcriptome assembly using the StringTie merge operation. The GffCompare program (Pertea and Pertea, 2020) was used to compare merged transcripts and reference annotations. The LncRNAs were identified using the following steps: (1) transcripts were retained as the potential lncRNAs with the class codes of “u,” “i,” “x,” and “o”; (2) transcripts were retained with a length of >200 bp; (3) the coding potential of each transcript was estimated using the coding potential calculator (CPC) (Kang et al., 2017), coding-non-coding index (CNCI, version 2) (Sun et al., 2013), and PfamScan (Mistry et al., 2013). Transcripts were identified as lncRNAs if they were marked as “non-coding” by both CPC and CNCI as well as failed to map to the Pfam database. The expression level of each lncRNA was normalized to FPKM.

We identified putative target genes of lncRNAs by integrating differentially expressed lncRNAs (DElncRNAs) and differentially expressed genes (DEGs) in pairwise comparison among stages, i.e., ST1 vs. ST2, ST2 vs. ST3-1, and ST3-1 vs. ST3-2. According to the cis-regulation mechanism, lncRNAs can regulate their upstream and downstream genes. Thus, DEGs located in 10 kb upstream and downstream of DElncRNAs were treated as candidate target genes of DElncRNAs. In contrast, DEGs trans-regulation by DElncRNAs were predicted using lncTar (Li et al., 2015). The regulatory networks of lncRNA-mRNA were visualized using Cytoscape (Saito et al., 2012).

Identification of miRNAs and miRNA Target Prediction

First, we converted clean reads from sRNA-seq into the required FASTA format of miRDeep-P2 (Kuang et al., 2019). Then, the FASTA files were mapped to the ginkgo genome using Bowtie (Langmead, 2010). The mapped reads were searched against the Rfam database (Kalvari et al., 2018; Kalvari et al., 2021) to remove rRNAs, tRNAs, snoRNAs, and snRNAs. The remaining reads were compared with conserved plant miRNAs in miRBase² (Release 22.1) (Kozomara and Griffiths-Jones, 2014) using Bowtie with no more than one mismatch. The remaining reads were used to predict novel miRNAs using miRDeep-P2. The expression level of each miRNA was normalized to reads per million mapped reads (RPM). The targets of miRNAs were predicted using psRNATarget (Dai et al., 2018).

Differential Expression Analysis and Function Enrichment Analysis

Differentially expressed genes, DElncRNAs, and differentially expressed miRNAs (DEmiRNAs) were identified by comparing two adjacent stages, i.e., ST1 vs. ST2, ST2 vs. ST3-1, and ST3-1 vs. ST3-2, with |log₂FC| > 1 and adjusted p-value < 0.05 using the R package DESeq2 (Love et al., 2014). To indicate the expression patterns of DEGs and DElncRNAs, we grouped DEGs and DElncRNAs into different clusters using the R package Mfuzz (Kumar and Futschik, 2007).

The genome of ginkgo was functionally annotated using eggNOG (Huerta-Cepas et al., 2017). Gene ontology (GO) annotation was analyzed using the R package clusterProfiler (Yu et al., 2012) with a p-value of 0.05.

Construction of ceRNA Regulatory Network

We constructed a lncRNA-miRNA-mRNA network based on the ceRNA hypothesis by integrating all DElncRNAs, DEmiRNAs, and DEGs in the ST2 vs. ST3-1 comparison as follows: (1) miRNA-lncRNA and miRNA-mRNA pairs were predicted using PsRNATarget with the default settings; (2) lncRNA-miRNA-mRNA pairs were integrated based on lncRNAs, and mRNAs were targeted by a common miRNA; and (3) the correlations of miRNA-lncRNA pairs and miRNA-mRNA pairs were evaluated based on the Spearman’s rank correlation coefficient (SCC) with ρ < -0.6. The Pearson correlation coefficient (PCC) was used to evaluate the strength of correlations between lncRNA-mRNA pairs, and pairs with r > 0.6 were selected. The lncRNA-miRNA-mRNA regulatory networks were visualized using Cytoscape.

Results

Morphology of Female Strobilus Buds During Early Differentiation

The differentiation process of ginkgo FSBs lasts nearly one year, which begins in May (for example, May 2020) and ends up in April of the following year (April 2021). This study focused on its early differentiation. The first stage of FSB differentiation was the vegetative stage when 3–7 leaf primordia curled inward in the buds before the end of May (ST1, Figure 1A). At the initial differentiation stage, an ovulate stalk primordium initiated at the axils of a leaf primordium or a bud bract, which was straighter than leaf primordia and in the form of a three-faced pyramid (ST2, Figure 1B). The subsequent stage is the exuberant differentiation stage (ST3), which was divided into two substages due to the long-lasting period from early July to early December. The first substage was characterized by the appearance and enlargement of multiple ovulate stalk primordia, and the ovulate stalk primordia were frontally laminar shape in appearance (ST3-1, Figure 1C). The second substage was featured by a longitudinal furrow appearing on the distal abaxial side of ovulated stalk primordium (ST3-2, Figure 1D).

FIGURE 1

Figure 1. Morphology of the female strobilus buds of ginkgo at three developmental stages was observed using scanning electron microscopy (SEM). (A) At the vegetative stage (ST1), 3–7 leaf primordia curl inward in the buds but no ovule stalk appeared. (B) At the initial differentiation stage (ST2), an ovulate stalk primordium is initiated at the axils of a leaf primordium. (C) At the early exuberant differentiation stage (ST3-1), multiple ovulate stalk primordia appear. With the enlargement of ovulate stalk primordia, some of them take on frontally laminar shapes (D). At the late exuberant differentiation stage (ST3-2), a longitudinal furrow becomes evident on the distal abaxial side of some ovulate stalk primordia. L, leaf primordium; OS, ovulate stalk primordium.

Expression Landscape of Genes During Early Differentiation of Female Strobilus Buds

Mapping clean reads from 12 strand-specific libraries to the ginkgo genome resulted in a 94.19–95.93% overall alignment rate. In total, 32,775 genes were expressed in at least one of the 12 samples. Furthermore, by comparing two adjacent stages (|log₂FC| > 1 and adjusted p-value < 0.05), 8,503 DEGs were identified, including 169 DEGs between ST1 and ST2, 6,515 DEGs between ST2 and ST3-1, and 1,819 DEGs between ST3-1 and ST3-2, respectively (Figure 2A and Supplementary Table 1). The greatest difference in gene expression indicated that associated genes changed dramatically from the initial differentiation stage to the exuberant differentiation stage (Figure 2A).

FIGURE 2

Figure 2. Characteristics of genes during the early differentiation of female strobilus buds in ginkgo. (A) The UpsetR plot shows the differentially expressed genes (DEGs) (|log₂FC| > 1; adjusted p-value < 0.05) between the vegetative (ST1) and the initial differentiation (ST2) stages, the initial differentiation (ST2), and the early exuberant differentiation (ST3-1) stages, the early exuberant differentiation (ST3-1) and the late exuberant differentiation (ST3-2) stages. (B) The fuzzy C-means clustering identified four distinct temporal patterns of DEGs. The x-axis represents the differentiation stages of ginkgo female strobilus buds, while the y-axis represents the normalized expression changes in each stage.

The fuzzy C-means clustering identified four distinct temporal patterns of DEGs (Figure 2B). Clusters 2 and 4 represented downregulated and upregulated genes, respectively, while clusters 1 and 3 suggested that genes showed bimodal expression patterns. Interestingly, Gb_19790 (GbFT) and Gb_13989 (LFY-like) showed the highest expression level in ST2 and were clustered into cluster 3. Moreover, Gb_16301 (AGAMOUS, AG) was also clustered into cluster 3, which plays a role in ovule identity and development (Pinyopich et al., 2003). Gb_32549 (CUP-SHAPED COTYLEDON 2, CUC2) in cluster 4 was differentially expressed between ST2 and ST3-1, and the highest expression level was found in ST3-2.

Gene ontology analysis indicated that gene functions differed among the four clusters (Supplementary Table 2). The upregulated genes in cluster 4 fell into the categories of meristem development, meristem initiation, and meristem structural organization, among others. In contrast, the downregulated genes in cluster 2 were enriched in 32 categories, including response to triterpenoid metabolic process, isopentenyl diphosphate biosynthetic process, mevalonate pathway, and response to jasmonic acid. Genes in cluster 1 exhibited the highest expression level at ST3-1, which were involved in the process of mitochondrial RNA modification and metabolic processes. The expression level of genes in cluster 3 reached the highest point at ST2, and GO enrichment analysis of these genes showed that the top GO terms were response to heat, secondary metabolic process, and response to water deprivation.

Expression Landscape of lncRNAs During Early Differentiation of Female Strobilus Buds

Besides protein-coding genes, 12,408 lncRNAs were identified based on the transcript length, the type of class code, and the coding potential from 12 strand-specific libraries. Up to 95.2% of lncRNAs were continuously expressed during the early differentiation of FSBs, while only a few lncRNAs were stage-specific (Figure 3A). The average length of lncRNAs was shorter than that of mRNAs (801.30 nt vs. 1298.07 nt, Figure 3B). Meanwhile, lncRNAs contained strikingly fewer exons than mRNAs, i.e., one exon for 73.8% of lncRNAs vs. ≥ 2 exons for 73.3% of mRNAs (Figure 3C). The GC content of lncRNAs was significantly lower than that of mRNAs (p-value < 2.22e^–16, Student’s t-test, Figure 3D). These results are consistent with previous studies on cucumber (He et al., 2020), poplar (Liu et al., 2019), and Chinese cabbage (Song et al., 2021).

FIGURE 3

Figure 3. Characteristics of long non-coding RNAs (lncRNAs) during the early differentiation of female strobilus buds in ginkgo. (A) The UpsetR plot shows the number of lncRNAs identified from each stage, i.e., the vegetative (ST1), the initial differentiation (ST2), and the exuberant differentiation stages (ST3-1 and ST3-2). (B–D) Comparisons of lncRNAs and mRNAs with respect to the transcript length, exon number, and GC content.

Pairwise comparisons between the three developmental stages suggested that 12, 4,075, and 160 lncRNAs were differentially expressed for ST1 vs. ST2, ST2 vs. ST3-1, and ST3-1 vs. ST3-2 comparisons, respectively (Supplementary Table 3). The fuzzy C-means algorithm also clustered all of the DElncRNAs into four clusters (Supplementary Figure 1). These four clusters showed a similar pattern to that of DEGs, implying that DElncRNAs may play an important role in the differentiation process of FSBs.

Long non-coding RNAs often mediate transcriptional regulation via cis- or trans-regulation (Rinn and Chang, 2012). As for cis-regulated analysis, 380 lncRNA-mRNA pairs, including 308 DEGs located 10 kb upstream and downstream of 351 DElncRNAs, composed the colocation network of the ST2 vs. ST3-1 comparison (Supplementary Figure 2). GO analysis showed that these DEGs were related to photosynthesis (Supplementary Table 4). In addition, 11 lncRNA-mRNA pairs, including 11 DElncRNAs and nine DEGs, formed the colocation network of ST3-1 vs. ST3-2 comparison (Supplementary Figure 3). Genes in this colocation network were significantly enriched in six GO categories, such as organic hydroxy compound biosynthetic process, saponin metabolic process, and glycoside biosynthetic process (Supplementary Table 4). In contrast, lncRNAs work in trans-regulation when they modulate genes across multiple chromosomes (Kornienko et al., 2013). We found that MSTRG.3148.1 modulated five DEGs in the ST2 vs. ST3-1 comparison by trans-regulation (Figure 4A), and MSTRG.39464.1 trans-regulated 21 DEGs, which composed a co-expression network in the ST3-1 vs. ST3-2 comparison (Figure 4B). However, no lncRNA was detected in regulating protein-coding genes through cis- and trans-regulation in the ST1 vs. ST2 comparison.

FIGURE 4

Figure 4. The lncRNA-mRNA co-expression networks during the initiation and enlargement of ovulate stalk primordia in ginkgo. (A) The lncRNA-mRNA co-expression network trans-regulated by differentially expressed lncRNAs (DElncRNAs) between the initial differentiation (ST2) and the early exuberant differentiation stages (ST3-1). (B) The lncRNA-mRNA co-expression network trans-regulated by DElncRNAs between the early (ST3-1) and the late (ST3-2) exuberant differentiation stages.

Expression Landscape of miRNAs During Early Differentiation of Female Strobilus Buds

The most mapped clean reads are 24 nt in length from 12 small RNA libraries (Figure 5A). We identified 398 conserved and 471 novel miRNAs totally. Most miRNAs were stage-specific, while only 239 miRNAs were expressed in all samples (Figure 5B), 57 DEmiRNAs (| log₂FC| > 1 and adjusted p-value < 0.05), including 11, 48, and 5 DEmiRNAs in the ST1 vs. ST2, ST2 vs. ST3-1, and ST2 vs. ST3-1 comparisons, respectively (Supplementary Table 5). Such a difference supported a remarkable transition from the initial differentiation to the early exuberant differentiation stage, in which miRNAs might play a role.

FIGURE 5

Figure 5. Characteristics of microRNAs (miRNAs) during the early differentiation of female strobilus buds in ginkgo. (A) Length distribution of the mappable sRNA reads from 12 libraries. (B) The UpsetR plot shows the number of miRNAs identified from each stage, i.e., the vegetative (ST1), the initial differentiation (ST2), and the exuberant differentiation stages (ST3-1 and ST3-2).

PsRNATarget analysis predicted that 144 DEGs were negatively targeted by 12 novel DEmiRNAs in the ST2 vs. ST3-1 comparison (Figure 6A) and eight DEGs were targeted by three novel DEmiRNAs in the ST3-1 vs. ST3-2 comparison (Figure 6B), while none of them were found in the ST1 vs. ST2 comparison. In the ST2 vs. ST3-1 comparison, three novel miRNAs, i.e., Gb_Novel_miR397, Gb_Novel_miR464, and Gb_Novel_miR189, were predicted to target 33, 18, and 18 genes, respectively. GO analysis showed that three genes (Gb_21271/Gb_21272/Gb_21273) targeted by Gb_Novel_miR189 and Gb_26565 targeted by Gb_Novel_miR397 were enriched in the GO term of positive regulation of cellular amide metabolic process. In the ST3-1 vs. ST3-2 comparison, Gb_Novel_miR331 might target six genes. However, target genes of known miRNAs were not found in this study.

FIGURE 6

Figure 6. The miRNA-mRNA regulatory networks during the initiation and enlargement of ovulate stalk primordia in ginkgo. (A) The target of genes of differentially expressed miRNAs (DEmiRNAs) between the initial differentiation (ST2) and the early exuberant differentiation stages (ST3-1). (B) The target of genes of differentially expressed miRNAs (DEmiRNAs) between the early (ST3-1) and the late (ST3-2) exuberant differentiation stages.

Long Non-coding RNA-Associated ceRNA Network Related to the Initiation and Enlargement of Ovulate Stalk Primordium

The ceRNA network based on DElncRNAs, DEmiRNAs, and DEGs in the ST3-1 vs. ST3-2 comparison predicted one lncRNA-miRNA-mRNA pair, i.e., MSTRG.2203.13-Gb_Novel_miR190- Gb_35828. In addition, 10,247 lncRNA-miRNA-mRNA pairs, which contained 1,225 nodes (393 lncRNAs, 44 miRNAs, and 788 mRNAs) and 1,361 edges (miRNA-target pairs), were found in the ST2 vs. ST3-1 comparison (Figure 7 and Supplementary Table 6). A total of three novel miRNAs, i.e., Gb_Novel_miR369, Gb_Novel_miR096, and Gb_Novel_miR397, were involved in most ceRNA pairs, suggesting that they may act as core regulators in the initiation and enlargement of ovulate stalk primordium in ginkgo. Of note, 13 lncRNAs were found to act as ceRNAs to compete for binding Gb_miR166b-5p and regulated the expression of 22 genes. The ceRNA pairs regulated by Gb_miR535-5p included three lncRNAs and four mRNAs. GO analysis showed that 23 genes that responded to water were regulated by 16 DEmiRNAs and 200 DElncRNAs (Figure 8). These results demonstrated the potential function of the ceRNA network during the initiation and enlargement of ovulate stalk primordium in ginkgo.

FIGURE 7

Figure 7. The long non-coding RNA-associated competing endogenous RNA (ceRNA) network between the early (ST3-1) and the late (ST3-2) exuberant differentiation stages in ginkgo.

FIGURE 8

Figure 8. A sub-competing endogenous RNA (ceRNA) network is involved in GO term of response to water between the early (ST3-1) and the late (ST3-2) exuberant differentiation stages in ginkgo.

Discussion

Protein-Coding Genes Are Dominant When the Differentiation of Female Strobilus Buds Starts

Two “flowering-time” genes, i.e., FT and SOC1, act as integrators in the regulatory network for the timing of floral transition in Arabidopsis. GbFT was verified to promote vegetative-to-reproductive induction (Wang et al., 2019), which is also supported in this study. The expression level of Gb_19790 (GbFT) was the highest at the initial differentiation stage and decreased along with the differentiation of strobilus buds. In contrast, SOC1, a member of the MIKC-type MADS-box genes, was named GbMADS6 (Gb_01884) and is believed to promote the transition of the vegetative shoot meristem into the reproductive meristem in ginkgo (Feng et al., 2021). In this study, the trend of its expression is the same as that of Feng et al. (2021), though GbMADS6 was not differentially expressed between the vegetative (early May) and the initial differentiation stages (end of May). The function of GbMADS6 needs to be further studied in future works.

LEAFY patterns the reproductive meristem by regulating APETALA3 (AP3, B function) and AGAMOUS (AG, C function) genes in angiosperms (Moyroud et al., 2011) and is conserved in gymnosperms and angiosperms (Moyroud et al., 2010). The two paralogous LFY-like genes in ginkgo, i.e., Gb_05859 (GinLfy) and Gb_13989 (GinNdly), showed different expression patterns that the expression level of GinNdly reached the highest at the initial differentiation stage, whereas GinLfy did not express differently in every pairwise comparison. We proposed that GinNdly was related to the development of FBSs in ginkgo, which is consistent with a previous study (Guo et al., 2005). The expression of WeILFY and/or WeINDLY always precedes that of WeIAP3/PI, and there is strong WeIGA expression in the male tissues of Welwitschia (Moyroud et al., 2017). Here, the expression pattern of Gb_16301 (AG) was the same as that of GinNdly (Figure 2B). Therefore, we proposed that controlling the AG-like MADS-box gene by NDLY homologous predated the origin of the female reproductive tissues in ginkgo.

Besides protein-coding genes, we identified 12,120 and 12,092 lncRNAs in the vegetative and the initial differentiation stages, respectively (Figure 3A). COOLAIR, a cold-induced FLC antisense transcript, played a role in early cold-induced silencing of FLC transcription (Swiezewski et al., 2009). However, we failed to predict COOLAIR, which may result from the FLC subfamily in MIKC-type genes having been lost in ginkgo (Chen et al., 2017; Liu et al., 2021). There were 12 lncRNAs and 11 miRNAs differently expressed, but their putative target genes were not predicted from our results. We proposed that protein-coding genes are dominant when the differentiation of ginkgo FSB starts.

Both Genes and Non-coding RNAs Are Involved in the Initiation and Enlargement of Ovulate Stalk Primordia

The greatest expression differences of genes, lncRNAs, and miRNAs were found between the initial differentiation and the early exuberant differentiation stages, i.e., 6,515 DEGs, 4,075 DElncRNAs, and 48 DEmiRNAs (Supplementary Tables 1, 3, 5), which implied that the largest changes occurred during the initiation and the enlargement of the ovulate stalk primordium. Of note, 1,819 DEGs, 160 DElncRNAs, and five DEmiRNAs were found between the early and late differentiation stages (Supplementary Tables 1, 3, 5). CUC proteins, the members of the NAC (NAM/ATAF/CUC) family, play important roles as boundary-defining factors during plant development in Arabidopsis (Ishida et al., 2000; Goncalves et al., 2015). Gb_32549 (CUC2) was upregulated at the early exuberant differentiation stage, which indicated that the expression of Gb_32549 can predict the boundary of the ovulated stalk primordium in ginkgo.

As for lncRNAs, cis- and trans-regulated analyses indicated that 362 and two DElncRNAs mediated the transcriptional regulation in the initiation and the enlargement of ovulate stalk primordia, respectively. The cis-regulating pairs were predominant in lncRNA-mRNA pairs here, which indicated that lncRNAs mainly regulated their target genes by cis-regulating in ginkgo (Supplementary Figures 2, 3). GO analysis showed that genes targeted by DElncRNAs had diverse functions (Supplementary Table 4). Three genes (Gb_09034, Gb_09035, and Gb_26289) cis-regulated by five lncRNAs were significantly enriched in the GO term of chlorophyll-binding. Gb_11209 and Gb_32759 were related to the biosynthetic and metabolic of the organic hydroxy compound and were cis-regulated by MSTRG.39998.1 and MSTRG.39391.1, respectively. In contrast, MSTRG.3148.1 and MSTRG.39464.1 functioned through trans-regulation and targeted five and 21 genes, respectively (Figure 4). Nine conserved miRNA families were related to flower development in plants. In this study, eight (miR156, miR159, miR160, miR166/165, miR167, miR169, miR172, and miR319) of them were identified (Luo et al., 2013), and 15 novel DEmiRNAs were predicted to target 152 DEGs (Figure 6). Gb_Novel_miR397 can target the greatest number of genes among DEmiRNAs, suggesting its essential roles in the initiation and enlargement of ovulate stalk primordia.

LncRNAs Act as ceRNAs Participating in the Initiation and Enlargement of Ovulate Stalk Primordia

A lncRNA-associated ceRNA network comprised of 10,248 lncRNA-miRNA-mRNA pairs was determined to be associated with the initiation and enlargement of ovulate stalk primordia (Figure 7 and Supplementary Table 6). Some genes and miRNAs were proven to regulate the development of flower buds and ovules previously. ARGONAUTE10 (AGO10) functions as a decoy for miR166/165 to regulate HD-ZIP III genes in SAM development and maintenance in Arabidopsis (Zhu et al., 2011). However, none of the genes targeted by Gb_miR166b-5p belongs to the HD-ZIP III transcription factors in this study. The function of Gb_miR166b-5p remains to be determined. MiR535 negatively regulated rice immune responses and cold tolerance in rice (Sun et al., 2020; Zhang et al., 2022); four genes (Gb_07340, Gb_10299, Gb_18723, Gb_30145) and three lncRNAs (MSTRG.14303.4, MSTRG.14303.9, MSTRG.23460.1) were predicted by Gb_miR535-5p. Water is confirmed to be related to the initiation and enlargement of ovule primordia in angiosperms (Clepet et al., 2021). Here, a sub-ceRNA network related to water consists of 200 lncRNAs, 16 miRNAs, and 23 genes. The ceRNA network may add more layers of complexity in the regulatory network of initiation and enlargement of ovulate stalk primordia.

Conclusion

As the first attempt, this study performed a genome-wide analysis to determine the roles of genes, lncRNAs, and miRNAs during the early differentiation of ginkgo FSBs. The temporal dynamic expression landscapes of genes and ncRNAs suggest that (1) the functions of some genes are conserved in both gymnosperms and angiosperms, such as GbFT (FT), GinNdly (LFY-like), Gb_16301(AG), and Gb_32549 (CUC2); (2) apart from genes, lncRNAs and miRNAs are closely associated with the dramatic morphological changes during the initiation and enlargement of ovulate stalk primordia; and (3) lncRNAs play vital roles as ceRNAs during the differentiation of flower buds. To conclude, this study provides key candidate transcripts related to the differentiation of strobilus buds for further wet-bench experiments in gymnosperms.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository and accession number can be found below: GinkgoDB with the accession number T0001 (https://ginkgo.zju.edu.cn/project/T0001/).

Author Contributions

Y-PZ conceived the ideas and supervised this project. P-PB, YS, and H-YL collected the samples and performed the experiments. P-PB conducted data analysis. P-PB and H-YL wrote the manuscript with the contributions from all co-authors. All authors read and approved the final manuscript.

Funding

This study was supported by the National Natural Science Foundation of China (Nos. 31870190 and 32071484).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors are grateful to Wen-Hao Li, Chen-Feng Lin, and Xin-Qi Liu for their help in collecting samples. They also thank Xin Zhang for constructive discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.863330/full#supplementary-material

Footnotes

1. ^ http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
2. ^ http://www.mirbase.org/

References

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Blumel, M., Dally, N., and Jung, C. (2015). Flowering time regulation in crops — what did we learn from Arabidopsis? Curr. Opin. Biotechnol. 32, 121–129. doi: 10.1016/j.copbio.2014.11.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, F., Zhang, X., Liu, X., and Zhang, L. (2017). Evolutionary analysis of MIKC(c)-type MADS-box genes in gymnosperms and angiosperms. Front. Plant Sci. 8:895. doi: 10.3389/fpls.2017.00895

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Chen, Y., Shi, C., Huang, Z., Zhang, Y., Li, S., et al. (2018). SOAPnuke: a mapreduce acceleration-supported software for integrated quality control and preprocessing of high-throughput sequencing data. Gigascience 7, 1–6. doi: 10.1093/gigascience/gix120

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, L. H., Yoon, J., and An, G. (2017). The control of flowering time by environmental factors. Plant J. 90, 708–719. doi: 10.1111/tpj.13461

PubMed Abstract | CrossRef Full Text | Google Scholar

Chuck, G., Cigan, A. M., Saeteurn, K., and Hake, S. (2007). The heterochronic maize mutant Corngrass1 results from overexpression of a tandem microRNA. Nat. Genet. 39, 544–549. doi: 10.1038/ng2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Clepet, C., Devani, R. S., Boumlik, R., Hao, Y. W., Morin, H., Marcel, F., et al. (2021). The miR166-SlHB15A regulatory module controls ovule development and parthenocarpic fruit set under adverse temperatures in tomato. Mol. Plant 14, 1185–1198. doi: 10.1016/j.molp.2021.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Dai, X. B., Zhuang, Z. H., and Zhao, P. X. C. (2018). psRNATarget: a plant small RNA target analysis server (2017 release). Nucleic. Acids Res. 46, W49–W54. doi: 10.1093/nar/gky316

PubMed Abstract | CrossRef Full Text | Google Scholar

Diggle, P. K., Di Stilio, V. S., Gschwend, A. R., Golenberg, E. M., Moore, R. C., Russell, J. R. W., et al. (2011). Multiple developmental processes underlie sex differentiation in angiosperms. Trends Genet. 27, 368–376. doi: 10.1016/j.tig.2011.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Dornelas, M. C., and Rodriguez, A. P. M. (2005). A FLORICAULA/LEAFY gene homolog is preferentially expressed in developing female cones of the tropical pine Pinus caribaea var. caribaea. Genet. Mol. Biol. 28, 299–307. doi: 10.1590/S1415-47572005000200021

CrossRef Full Text | Google Scholar

Douglas, A. W., Stevenson, D. W., and Little, D. P. (2007). Ovule development in Ginkgo biloba L., with emphasis on the collar and nucellus. Int. J. Plant Sci. 168, 1207–1236. doi: 10.1086/521693

CrossRef Full Text | Google Scholar

Feng, Z., Yang, T. T., Li, M., Dong, J. J., Wang, G. B., Wang, Q. Y., et al. (2021). Identification and cloning of GbMADS6, a SOC1 homolog gene involved in floral development in Ginkgo biloba. J. Plant Biochem. Biotechnol. 30, 554–563. doi: 10.1007/s13562-021-00646-4

CrossRef Full Text | Google Scholar

Fornara, F., de Montaigu, A., and Coupland, G. (2010). SnapShot: control of flowering in Arabidopsis. Cell 141, 550–550. doi: 10.1016/j.cell.2010.04.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Frohlich, M. W., and Parker, D. S. (2000). The mostly male theory of flower evolutionary origins: from genes to fossils. Syst. Bot. 25, 155–170. doi: 10.2307/2666635

CrossRef Full Text | Google Scholar

Fu, D. Z., and Yang, Q. E. (1993). A new morohological interpretation of the female reproductive organs in Ginkgo biloba., with a phyligenetic considertation on gymnosperms. Acta Phytotaxon. Sin. 31, 309–317.

Google Scholar

Goncalves, B., Hasson, A., Belcram, K., Cortizo, M., Morin, H., Nikovics, K., et al. (2015). A conserved role for CUP-SHAPED COTYLEDON genes during ovule development. Plant J. 83, 732–742. doi: 10.1111/tpj.12923

PubMed Abstract | CrossRef Full Text | Google Scholar

Guan, R., Zhao, Y., Zhang, H., Fan, G., Liu, X., Zhou, W., et al. (2019). Updated draft genome assembly of Ginkgo biloba. Gigascience Database. doi: 10.5524/100613

CrossRef Full Text | Google Scholar

Guan, R., Zhao, Y. P., Zhang, H., Fan, G. Y., Liu, X., Zhou, W. B., et al. (2016). Draft genome of the living fossil Ginkgo biloba. Gigascience 5:49. doi: 10.1186/s13742-016-0154-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, C.-L., Chen, L.-G., He, X.-H., Dai, Z., and Yuan, H.-Y. (2005). Expressions of LEAFY homologous genes in different organs and stages of Ginkgo biloba. Yichuan 27, 241–244.

PubMed Abstract | Google Scholar

He, C., Cai, Y., Li, M., Dong, J., Liu, W., Wang, J., et al. (2018). Screening and expression analysis of related genes based on transcriptome sequencing of ginkgo flower buds at three differentiation stages. Acta Hortic. Sin. 45, 1479–1490.

Google Scholar

He, X., Guo, S., Wang, Y., Wang, L., Shu, S., and Sun, J. (2020). Systematic identification and analysis of heat-stress-responsive lncRNAs, circRNAs and miRNAs with associated co-expression and ceRNA networks in cucumber (Cucumis sativus L.). Physiol. Plant 168, 736–754. doi: 10.1111/ppl.12997

PubMed Abstract | CrossRef Full Text | Google Scholar

Huerta-Cepas, J., Forslund, K., Coelho, L. P., Szklarczyk, D., Jensen, L. J., von Mering, C., et al. (2017). Fast genome-wide functional annotation through orthology assignment by eggNOG-Mapper. Mol. Biol. Evol. 34, 2115–2122. doi: 10.1093/molbev/msx148

PubMed Abstract | CrossRef Full Text | Google Scholar

Ishida, T., Aida, M., Takada, S., and Tasaka, M. (2000). Involvement of CUP-SHAPED COTYLEDON genes in gynoecium and ovule development in Arabidopsis thaliana. Plant Cell Physiol. 41, 60–67. doi: 10.1093/pcp/41.1.60

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalvari, I., Nawrocki, E. P., Argasinska, J., Quinones-Olvera, N., Finn, R. D., Bateman, A., et al. (2018). Non-coding RNA analysis using the Rfam database. Curr. Protoc. Bioinform. 62:e51. doi: 10.1002/cpbi.51

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalvari, I., Nawrocki, E. P., Ontiveros-Palacios, N., Argasinska, J., Lamkiewicz, K., Marz, M., et al. (2021). Rfam 14: expanded coverage of metagenomic, viral and microRNA families. Nucleic. Acids Res. 49, D192–D200. doi: 10.1093/nar/gkaa1047

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, Y. J., Yang, D. C., Kong, L., Hou, M., Meng, Y. Q., Wei, L. P., et al. (2017). CPC2: a fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res. 45, W12–W16. doi: 10.1093/nar/gkx428

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D., Paggi, J. M., Park, C., Bennett, C., and Salzberg, S. L. (2019). Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 37, 907–915. doi: 10.1038/s41587-019-0201-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Kornienko, A. E., Guenzl, P. M., Barlow, D. P., and Pauler, F. M. (2013). Gene regulation by the act of long non-coding RNA transcription. BMC Biol. 11:59. doi: 10.1186/1741-7007-11-59

PubMed Abstract | CrossRef Full Text | Google Scholar

Kozomara, A., and Griffiths-Jones, S. (2014). miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic. Acids Res 42, D68–D73. doi: 10.1093/nar/gkt1181

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuang, Z., Wang, Y., Li, L., and Yang, X. (2019). miRDeep-P2: accurate and fast analysis of the microRNA transcriptome in plants. Bioinformatics 35, 2521–2522. doi: 10.1093/bioinformatics/bty972

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, L., and Futschik, M. E. (2007). Mfuzz: a software package for soft clustering of microarray data. Bioinformation 2, 5–7. doi: 10.6026/97320630002005

PubMed Abstract | CrossRef Full Text | Google Scholar

Langmead, B. (2010). Aligning short sequencing reads with Bowtie. Curr. Protoc. Bioinformatics Chapter 11, Unit 11.7. doi: 10.1002/0471250953.bi1107s32

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Ma, W., Zeng, P., Wang, J., Geng, B., Yang, J., et al. (2015). LncTar: a tool for predicting the RNA targets of long noncoding RNAs. Brief. Bioinform. 16, 806–812. doi: 10.1093/bib/bbu048

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H., Wang, X., Wang, G., Cui, P., Wu, S., Ai, C., et al. (2021). The nearly complete genome of Ginkgo biloba illuminates gymnosperm evolution. Nat. Plants 7, 748–756. doi: 10.1038/s41477-021-00933-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, S., Wu, L., Qi, H. R., and Xu, M. (2019). LncRNA/circRNA-miRNA-mRNA networks regulate the development of root and shoot meristems of Populus. Ind. Crop. Prod. 133, 333–347. doi: 10.1016/j.indcrop.2019.03.048

CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, Y., Guo, Z., and Li, L. (2013). Evolutionary conservation of microRNA regulatory programs in plant flower development. Dev. Biol. 380, 133–144. doi: 10.1016/j.ydbio.2013.05.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, D., Ye, J. B., and Xu, F. (2019). Advances of the flowering genes of gymnosperms. Not. Bot. Horti Agrobot. Cluj-Na. 47, 1–9. doi: 10.15835/nbha47111343

CrossRef Full Text | Google Scholar

Mistry, J., Finn, R. D., Eddy, S. R., Bateman, A., and Punta, M. (2013). Challenges in homology search: HMMER3 and convergent evolution of coiled-coil regions. Nucleic Acids Res. 41:e121. doi: 10.1093/nar/gkt263

PubMed Abstract | CrossRef Full Text | Google Scholar

Moon, J., Suh, S. S., Lee, H., Choi, K. R., Hong, C. B., Paek, N. C., et al. (2003). The SOC1 MADS-box gene integrates vernalization and gibberellin signals for flowering in Arabidopsis. Plant J. 35, 613–623. doi: 10.1046/j.1365-313X.2003.01833.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Moyroud, E., Kusters, E., Monniaux, M., Koes, R., and Parcy, F. (2010). LEAFY blossoms. Trends Plant Sci. 15, 346–352. doi: 10.1016/j.tplants.2010.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Moyroud, E., Minguet, E. G., Ott, F., Yant, L., Pose, D., Monniaux, M., et al. (2011). Prediction of regulatory interactions from genome sequences using a biophysical model for the Arabidopsis LEAFY transcription factor. Plant Cell 23, 1293–1306. doi: 10.1105/tpc.111.083329

PubMed Abstract | CrossRef Full Text | Google Scholar

Moyroud, E., Monniaux, M., Thevenon, E., Dumas, R., Scutt, C. P., Frohlich, M. W., et al. (2017). A link between LEAFY and B-gene homologues in welwitschia mirabilis sheds light on ancestral mechanisms prefiguring floral development. New Phytol. 216, 469–481. doi: 10.1111/nph.14483

PubMed Abstract | CrossRef Full Text | Google Scholar

Navarro, C., Abelenda, J. A., Cruz-Oro, E., Cuellar, C. A., Tamaki, S., Silva, J., et al. (2011). Control of flowering and storage organ formation in potato by FLOWERING LOCUS T. Nature 478, 119–122. doi: 10.1038/nature10431

PubMed Abstract | CrossRef Full Text | Google Scholar

Pertea, G., and Pertea, M. (2020). GFF Utilities: GffRead and GffCompare. F1000Res 9:304. doi: 10.12688/f1000research.23297.2

PubMed Abstract | CrossRef Full Text | Google Scholar

Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T. C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295. doi: 10.1038/nbt.3122

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinyopich, A., Ditta, G. S., Savidge, B., Liljegren, S. J., Baumann, E., Wisman, E., et al. (2003). Assessing the redundancy of MADS-box genes during carpel and ovule development. Nature 424, 85–88. doi: 10.1038/nature01741

PubMed Abstract | CrossRef Full Text | Google Scholar

Rinn, J. L., and Chang, H. Y. (2012). Genome regulation by long noncoding RNAs. Annu Rev. Biochem. 81, 145–166. doi: 10.1146/annurev-biochem-051410-092902

PubMed Abstract | CrossRef Full Text | Google Scholar

Saito, R., Smoot, M. E., Ono, K., Ruscheinski, J., Wang, P. L., Lotia, S., et al. (2012). A travel guide to Cytoscape plugins. Nat. Methods 9, 1069–1076. doi: 10.1038/nmeth.2212

PubMed Abstract | CrossRef Full Text | Google Scholar

Salmena, L., Poliseno, L., Tay, Y., Kats, L., and Pandolfi, P. P. (2011). A ceRNA hypothesis: the rosetta stone of a hidden RNA language? Cell 146, 353–358. doi: 10.1016/j.cell.2011.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, J. K., Fan, W. G., and Wen, X. P. (1998). Female flower bud diferentiation of Ginkgo biloba. Acta Hortic. Sin. 25, 33–36.

Google Scholar

Shindo, S., Sakakibara, K., Sano, R., Ueda, K., and Hasebe, M. (2001). Characterization of a FLORICAULA/LEAFY homologue of Gnetum parvifolium and its implications for the evolution of reproductive organs in seed plants. Int. J. Plant Sci. 162, 1199–1209. doi: 10.1086/323417

CrossRef Full Text | Google Scholar

Song, X., Hu, J., Wu, T., Yang, Q., Feng, X., Lin, H., et al. (2021). Comparative analysis of long noncoding RNAs in angiosperms and characterization of long noncoding RNAs in response to heat stress in Chinese cabbage. Hortic. Res. 8:48. doi: 10.1038/s41438-021-00484-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, L., Luo, H., Bu, D., Zhao, G., Yu, K., Zhang, C., et al. (2013). Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic. Acids Res. 41, e166. doi: 10.1093/nar/gkt646

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, M. Z., Shen, Y., Yang, J. K., Cai, X. X., Li, H. Y., Zhu, Y. M., et al. (2020). miR535 negatively regulates cold tolerance in rice. Mol. Breed 40:14. doi: 10.1007/s11032-019-1094-0

CrossRef Full Text | Google Scholar

Swiezewski, S., Liu, F. Q., Magusin, A., and Dean, C. (2009). Cold-induced silencing by long antisense transcripts of an Arabidopsis Polycomb target. Nature 462, 799–802. doi: 10.1038/nature08618

PubMed Abstract | CrossRef Full Text | Google Scholar

Tandre, K., Albert, V. A., Sundas, A., and Engstrom, P. (1995). Conifer homologs to genes that control floral development in Angiosperms. Plant Mol.Biol. 27, 69–78. doi: 10.1007/Bf00019179

PubMed Abstract | CrossRef Full Text | Google Scholar

Theissen, G., Becker, A., Di Rosa, A., Kanno, A., Kim, J. T., Munster, T., et al. (2000). A short history of MADS-box genes in plants. Plant Mol. Biol. 42, 115–149. doi: 10.1023/A:1006332105728

CrossRef Full Text | Google Scholar

Theissen, G., and Saedler, H. (2001). Plant biology. Floral quartets. Nature 409, 469–471. doi: 10.1038/35054172

PubMed Abstract | CrossRef Full Text | Google Scholar

Walas, L., Mandryk, W., Thomas, P. A., Tyrala-Wierucka, Z., and Iszkulo, G. (2018). Sexual systems in gymnosperms: a review. Basic Appl. Ecol. 31, 1–9. doi: 10.1016/j.baae.2018.05.009

CrossRef Full Text | Google Scholar

Walden, A. R., Wang, D. Y., Walter, C., and Gardner, R. C. (1998). A large family of TM3 MADS-box cDNAs in Pinus radiata includes two members with deletions of the conserved K domain. Plant Sci. 138, 167–176. doi: 10.1016/S0168-9452(98)00160-5

CrossRef Full Text | Google Scholar

Wang, L. L., Yan, J. P., Zhou, X., Cheng, S. Y., Chen, Z. X., Song, Q. L., et al. (2019). GbFT, a FLOWERING LOCUS T homolog from Ginkgo biloba, promotes flowering in transgenic Arabidopsis. Sci. Hortic. 247, 205–215. doi: 10.1016/j.scienta.2018.12.020

CrossRef Full Text | Google Scholar

Weigel, D., Alvarez, J., Smyth, D. R., Yanofsky, M. F., and Meyerowitz, E. M. (1992). LEAFY controls floral meristem identity in Arabidopsis. Cell 69, 843–859. doi: 10.1016/0092-8674(92)90295-n

CrossRef Full Text | Google Scholar

Winter, K. U., Becker, A., Munster, T., Kim, J. T., Saedler, H., and Theissen, G. (1999). MADS-box genes reveal that gnetophytes are more closely related to conifers than to flowering plants. Proc. Natl. Acad. Sci. U.S.A. 96, 7342–7347. doi: 10.1073/pnas.96.13.7342

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, G., and Poethig, R. S. (2006). Temporal regulation of shoot development in Arabidopsis thaliana by miR156 and its target SPL3. Development 133, 3539–3547. doi: 10.1242/dev.02521

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Z. C., Yang, C. C., Wang, Z. Y., Yang, Z., Chen, D. Y., and Wu, Y. J. (2019). LncRNA expression profile and ceRNA analysis in tomato during flowering. PLoS One 14:e0210650. doi: 10.1371/journal.pone.0210650

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J. Y., Chen, L. G., Hu, X. Q., and He, X. H. (2002). GinNdly gene cloned from the male gametophyte of Ginkgo biloba. Chin. J. Cell Biol. 24, 189–191.

Google Scholar

Zhang, L. L., Huang, Y. Y., Zheng, Y. P., Liu, X. X., Zhou, S. X., Yang, X. M., et al. (2022). Osa-miR535 targets SQUAMOSA promoter binding protein-like 4 to regulate blast disease resistance in rice. Plant J. 8, doi: 10.1111/tpj.15663

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, T., Wang, J., and Zhou, C. (2015). The role of miR156 in developmental transitions in Nicotiana tabacum. Sci. China Life Sci. 58, 253–260. doi: 10.1007/s11427-015-4808-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y. C., Liao, J. Y., Li, Z. Y., Yu, Y., Zhang, J. P., Li, Q. F., et al. (2014). Genome-wide screening and functional analysis identify a large number of long noncoding RNAs involved in the sexual reproduction of rice. Genome Biol. 15:512. doi: 10.1186/s13059-014-0512-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, H., Hu, F., Wang, R., Zhou, X., Sze, S. H., Liou, L. W., et al. (2011). Arabidopsis argonaute10 specifically sequesters miR166/165 to regulate shoot apical meristem development. Cell 145, 242–256. doi: 10.1016/j.cell.2011.03.024

PubMed Abstract | CrossRef Full Text | Google Scholar