Skip to main content

Identification of the specific long-noncoding RNAs involved in night-break mediated flowering retardation in Chenopodium quinoa

Abstract

Background

Night-break (NB) has been proven to repress flowering of short-day plants (SDPs). Long-noncoding RNAs (lncRNAs) play key roles in plant flowering. However, investigation of the relationship between lncRNAs and NB responses is still limited, especially in Chenopodium quinoa, an important short-day coarse cereal.

Results

In this study, we performed strand-specific RNA-seq of leaf samples collected from quinoa seedlings treated by SD and NB. A total of 4914 high-confidence lncRNAs were identified, out of which 91 lncRNAs showed specific responses to SD and NB. Based on the expression profiles, we identified 17 positive- and 7 negative-flowering lncRNAs. Co-expression network analysis indicated that 1653 mRNAs were the common targets of both types of flowering lncRNAs. By mapping these targets to the known flowering pathways in model plants, we found some pivotal flowering homologs, including 2 florigen encoding genes (FT (FLOWERING LOCUS T) and TSF (TWIN SISTER of FT) homologs), 3 circadian clock related genes (EARLY FLOWERING 3 (ELF3), LATE ELONGATED HYPOCOTYL (LHY) and ELONGATED HYPOCOTYL 5 (HY5) homologs), 2 photoreceptor genes (PHYTOCHROME A (PHYA) and CRYPTOCHROME1 (CRY1) homologs), 1 B-BOX type CONSTANS (CO) homolog and 1 RELATED TO ABI3/VP1 (RAV1) homolog, were specifically affected by NB and competed by the positive and negative-flowering lncRNAs. We speculated that these potential flowering lncRNAs may mediate quinoa NB responses by modifying the expression of the floral homologous genes.

Conclusions

Together, the findings in this study will deepen our understanding of the roles of lncRNAs in NB responses, and provide valuable information for functional characterization in future.

Peer Review reports

Background

Plants sense environmental changes to maximize flowering transition success. Accumulating evidences have revealed that six major factors, including temperature, photoperiod, gibberellin acids (GAs), age and autonomous pathways are implicated in flowering regulation [1]. The photoperiod pathway controls floral transition through a signal cascade involving circadian clock and florigen genes. Under inductive photoperiods when the expression phases of endogenous circadian clock regulators coincide with external light signal, the clock output genes are converged to activate the florigen-encoding genes, FLOWERING LOCUS T (FT) and its homologs [2,3,4]. Thereafter, FT proteins are transported to shoot apical meristems where interact with FLOWERING LOCUS T (FD) and 14–3-3 proteins to form florigen activation complex (FAC), which in turn induces the meristem identity gene APETALA1 and triggers flowering [5, 6]. Day length is a critical factor for photoperiodic flowering. In short-day plants (SDPs) such as rice (Oryza saltiva) and soybean (Glycine max), FT transcripts are more abundant when the night length longer than a certain threshold [7, 8]. By contrast, in long-day plants (LDPs) such as Arabidopsis (Arabidopsis thaliana) and wheat (Triticum aestivum), only when the day length exceeds a certain threshold, FT proteins are sufficient for flowering induction [9, 10]. Different photoperiods combined with short pulses of light during the night period (referred to as night break, NB) approaches are well-established and widely adopted to study the photoperiodic regulation of flowering in many plant species. NB causes different effects in flowering of SDPs and LDPs by changing the expression of FT. In rice, Heading date 3a (Hd3a), a homolog of Arabidopsis FT, is the principle mediator responsible for NB flowering retardation [11]. In soybean, down-regulation of FT2a and FT5a is the major cause for NB responses [12]. In contrast, multiple NBs accelerate heading in wheat plants grown under SD by inducing FT1 [13]. Genetic studies indicate that the transcriptional changes of FT in response to NB are predominantly mediated by phytochromes. NB responses are abolished in the phyB mutants of rice [14], phyA mutants of soybean [12] and the phyB and phyC mutants of wheat [13], indicating functional conservation and divergence of phytochromes across different plant species. Further studies indicate circadian clock genes are also involved in NB responses. The up-regulation of FT1 by NB in wheat is dependent on PHOTOPERIOD1 (PPD1), a homolog of PSEUDO RESPONSE REGULATOR 37 (PRR37) in rice, which is a core component of the circadian clock [13]. In poplar (Populus L.), NB induces FT2 to free the shoot growth cessation by suppressing the circadian clock gene LATE ELONGATED HYPOCOTYL 2 (LHY2) [15]. Thus, NB affects plant flowering by controlling different layers of regulators.

Long-noncoding RNAs (lncRNAs) refer to a class of RNA transcripts longer than 200 nt which lack discernable protein coding potential [16, 17]. Compared with mRNAs, lncRNAs have much lower expression levels and sequence conservation [16, 17]. Most lncRNAs harbor spatio-temporal specificity [16, 17]. lncRNAs could be generated from the intergenic, intronic, or coding regions in the sense and antisense directions of mRNA [16, 17]. lncRNAs control the expression of target genes in cis- or trans-acting way [18, 19]. lncRNAs exert their functions by serving as sponge for miRNAs, functioning as precursors of miRNAs in the cytoplasm [20], or serving as scaffolds of epigenetic regulators to modulate chromatin status in the nucleus [17, 21]. A few studies indicate that lncRNAs participate in different biological processes, such as flowering, fertility, phosphate metabolism, leaf patterning, nodule formation and phytohormone-related development. COLD INDUCED LONG ANTISENSE INTRAGENIC RNAs (COOLAIR) and COLD ASSISTED INTRONIC NONCODING RNA (COLDAIR), two different classes of lncRNAs transcribed from the site of FLOWERING LOCUS C (FLC) in the antisense direction and from the first intron of FLC in the sense direction, respectively, are involved in flowering regulation by repressing FLC via epigenetic modulation [22,23,24]. LD-SPECIFIC MALE-FERTILITY-ASSOCIATED RNA (LDMAR) controls the photoperiod-sensitive male sterility in rice under LD [25]. The rice lncRNA INDUCED BY PHOSPHATE STARVATION 1 (IPS1) functions as a target mimic of miR399 to up-regulate the phosphate metabolism gene PHOSPHATE2 (PHO2) when encountered with phosphate starvation [26]. TWISTED LEAF (TL), an antisense lncRNA of OsMYB60, suppresses OsMYB60 via epigenetic chromatin modification to regulate leaf morphology [27]. In Medicago truncatula, the lncRNA EARLY NODULIN 40 (ENOD40) ENOD40 interacts with RNA-BINDING PROTEIN 1 (MtRBP1) to change its location from nuclear speckles to cytoplasmic granules to promote nodule formation [28]. The lincRNA AUXIN REGULATED PROMOTER LOOP (APOLO) regulates the promotor chromatin loop of PINOID (PID) which regulates polar auxin transport [29]. Recently, genome-wide identification of lncRNAs has been performed in a few plants. Huang et al identified 14 lncRNAs co-expressed with 10 pollen-associated coding genes during pollen development in Brassica rapa (B. rapa) [30]. One hundred forty-two lncRNAs may participate in fruit ripening and the climacteric in Cucumis melo by regulating the targets involved in auxin signal transduction, ethylene, sucrose biosynthesis and signaling and the ABA signaling pathway [31]. One hundred eighty-five salt stress-related lncRNAs were obtained in duckweed (Spirodela polyrhiza) [32]. By co-expression analysis, the important lncRNA-mRNA pairs associated with contrasting drought stress responses in two rapeseed cultivators (B. napus) were obtained [33]. These studies have shed lights on the roles of lncRNAs. However, the roles of lncRNAs in flowering regulation remain largely unknown, especially in quinoa (Chenopodium quinoa), an important SD coarse cereal.

Quinoa has been recommended by FAO as “super food” because of the nutritious elements in the seeds [34]. As quinoa is a facultative SD plant [35, 36], many photoperiod-sensitive quinoa accessions may encounter with yield penalties when grown across different latitudes. However, for the time being the lncRNAs and corresponding regulatory mechanisms involved in quinoa photoperiodic flowering remain to be elucidated. To this end, in the present study, we investigated the effect of NB on quinoa flowering, and by performing strand-specific RNA sequencing (SS RNA-seq) of the plants grown under NB and SD, we identified the specific coding genes and lncRNAs that may be involved in flowering regulation. Hence, this study provides novel insights into the quinoa flowering molecular mechanisms.

Results

Effects of short-day and night-break on quinoa flowering time

To evaluate the effects of short-day (SD) and night-break (NB) on quinoa flowering, two quinoa cultivars (“HL1” and “HL2”) with different maturation periods were used. Seeds were sown in pots and grown under natural day conditions before photoperiodic treatment. Then, 14 d-old seedlings were transferred into growth chambers and subjected to SD and NB treatment, respectively. After 14 d photoperiodic entertainment in growth chambers (28 d-after-sowing, 28 DAS), the floral buds initiated in “HL1” and “HL2” under SD, whereas no visible floral buds found in both cultivars under NB. On 40 DAS, the floral buds size was obviously bigger under SD (Fig. 1a, c) than that under NB (Fig. 1b, d). Thus, we speculated that SD accelerates whereas NB represses floral buds development.

Fig. 1
figure1

Phenotypes of the floral buds in two quinoa cultivars under SD and NB. On 40 DAS, the floral buds were apparently visible in cultivars “HL1” (a) and “HL2” (c) under SD, whereas under NB, the floral buds were very tiny in cultivars “HL1” (b) and “HL2” (d)

Transcriptome sequencing and identification of lncRNAs in quinoa

To uncover the roles of lncRNAs in NB responses of quinoa, seedlings of “HL1” were grown in growth chamber and treated by SD and NB as described in section Methods. The leaf samples of “HL1” were collected on 14 DAS (SD1), 26 DAS (SD2) under SD, and 28 DAS under NB. Sampling was performed with three biological replicates for each time point. Then the samples were subjected to total RNA extraction and transcriptome analysis. Through SS RNA-seq on an Illumina Hiseq4000 platform, an average of more than 16.02 GB data of 9 cDNA libraries were obtained (Table 1). The Q30 value was more than 90.79% for each library (Table 1). Subsequently, the clean reads were mapped to the quinoa reference genome [37] using HISAT2 software [38], and the mapping rate was higher than 85.48% for each library (Table 1). Based on the mapping results, StringTie [39] was used to re-construct the transcripts and predict the novel genes. The assembled transcripts were annotated by Gffcompare program [40], and the unannotated transcripts were further used for long-noncoding RNA (lncRNA) prediction. Based on the criteria described in section Methods, a total of 4914 lncRNAs were identified. According to their genomic origins, 4075 intergenic lncRNAs (lincRNAs), 250 intronic lncRNA, 360 sense lncRNAs, and 229 antisense lncRNAs were identified (Fig. 2a), respectively. By performing BLAST on lncRNAs database, we found only 605 lncRNAs were annotated in the NONCODE database [41], and most (87%) of the quinoa lncRNAs were novel genes. Meanwhile, as shown in Fig. 2b, various types of lncRNAs were evenly distributed over the whole quinoa genome at the density of 3.5 per Mb.

Table 1 Summary of the clean reads obtained from 9 strand-specific cDNA libraries
Fig. 2
figure2

Characteristics of the short-day and night-break responsive lncRNAs identified from quinoa leaves. a Numbers of four types of lncRNAs. b Distribution of different-type lncRNAs on the quinoa chromosomes. Five circles from outside to inside represents quinoa chromosomes, sense lncRNAs (green), intergenic lncRNAs (red), intron lncRNAs (blue), antisense lncRNAs (gray). c Comparison of the expression levels of lncRNAs and mRNAs. d Open reading frame (ORF) length distributions of lncRNA and mRNA transcripts

The relative gene expression level of each mRNA and lncRNA was normalized to Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced (FPKM) value [42]. Based on the gene expression profiles, Spearman correlation analysis suggested that, three biological replicates were highly related with correlation efficiencies larger than 0.94 (R2 > 0.94), and were readily clustered into the same clade (Additional file 1: Fig. S1). Besides, SD1 and NB samples were clustered as neighboring clades, distant from SD2 clade (Fig. S1), indicating some transcriptional changes from SD1 to SD2 were recovered by NB. Statistics of mRNA and lncRNA expression profiles showed that lncRNAs were expressed at much lower levels compared with mRNAs (Fig. 2c). The open reading frame (ORF) lengths comparison results showed that, lncRNAs predominately harbored 100 bp-long ORFs, significantly shorter than that of mRNAs (Fig. 2d). Based on a criterion of average FPKM of three replicates> 0.1, we found 1698 lncRNAs were commonly expressed over different time points, and 361, 569 and 415 lncRNAs were specifically expressed at SD1, SD2, and NB, respectively (Fig. 3a). Meanwhile, based on the threshold of average FPKM of three replicates> 1, 14,056 mRNAs were detected across the whole time points (Fig. 3b). More genes (762) were specifically expressed at SD2 compared with the expressed genes at SD1 (173) and NB (605) (Fig. 3b). Coincidently, more lncRNAs and mRNAs were specifically expressed at SD2 than at NB, indicating SD may promote the expression of flowering genes but NB may repress this effect.

Fig. 3
figure3

Venn diagrams showing the numbers of commonly and specifically expressed lncRNAs and mRNAs in different groups. a Numbers of commonly and specifically expressed lncRNAs in SD1, SD2 and NB based on a cutoff of FPKM> 0.1. b Numbers of commonly and specifically expressed mRNAs in SD1, SD2 and NB based on a cutoff of FPKM> 1. The FPKM value was denoted as the average value of three replicates

Differential expression analysis of lncRNAs and mRNAs

DESeq2 [43] was employed to explore the differentially expressed lncRNAs (DE lncRNAs) and mRNAs (DEGs) between different groups with thresholds of FDR < 0.05 and |log2(fold change)| ≥ 1. In total, 91 DE lncRNAs and 1401 DEGs were obtained in comparisons of SD1_vs_SD2 and SD2_vs_NB (Table 2, Additional file 2: Table S1, Additional file 3: Table S2). Sixty-four lncRNAs were differentially expressed in SD1_vs_SD2, out of which 40 and 24 lncRNAs were up- and down-regulated in SD1_vs_SD2 (Table 2), respectively. After transferring from SD2 to NB, the expression levels of 20 and 31 lncRNAs were increased and decreased, respectively (Table 2). DEGs analysis showed that 693 and 309 genes were up- and down-regulated in SD1_vs_SD2, and 315 and 379 were up- and down-regulated in SD2_vs_NB, respectively (Table 2). These results indicated that many lncRNA and mRNA genes tended to be activated by SD whereas repressed by NB.

Table 2 Brief review of the differentially expressed lncRNAs and mRNA between different groups

Functional annotation of the differentially expressed lncRNAs in response to SD and NB

To investigate the potential roles of DE lncRNAs, we performed functional annotation of the DE lncRNAs targets. The cis targets were selected from the mRNAs located at less than 100kbp downstream and upstream of each lncRNA. The trans targets were determined with cutoffs of Pearson’s correlation coefficient > 0.9 or < − 0.9 and P value < 0.01 between the expression profiles of mRNAs and lncRNAs. Co-expression and co-location analysis demonstrated that 5917 targets were regulated by the up-regulated lncRNAs in SD1_vs_SD2, including 29,401 lncRNA-mRNA pairs (Additional file 4: Table S3). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis indicated that these targets were abundant in pathways of “Ribosome biogenesis in eukaryotes” and “Photosynthesis-antenna proteins” (Fig. 4a). 5100 targets and 9340 lncRNA-mRNA pairs were predicted to be related with the down-regulated lncRNAs in SD1_vs_SD2 (Additional file 5: Table S4). “Photosynthesis-antenna proteins” was the only significantly enriched KEGG pathway for these target genes (Fig. 4b). In the network of the up-regulated lncRNAs in SD2_vs_NB, 4736 targets and 8420 lncRNA-mRNA pairs were contained (Additional file 6: Table S5). KEGG pathways including “Photosynthesis-antenna proteins” and “Glyoxylate and dicarboxylate metabolism” were highly related with these targets (Fig. 4c). 5929 targets and 20,438 lncRNA-mRNA pairs were involved in the network of the down-regulated lncRNAs in SD2_vs_NB (Additional file 7: Table S6). “Ribosome biogenesis in eukaryotes” and “Photosynthesis-antenna proteins” were the most significantly enriched KEGG pathways for these targets (Fig. 4d).

Fig. 4
figure4

Functional annotation of the differentially expressed lncRNAs by Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. KEGG enrichment analysis of the up-regulated lncRNA-targets (a) and down-regulated lncRNA-targets (b) in comparison of SD1_vs_SD2. KEGG enrichment analysis of the up-regulated lncRNA-targets (c) and down-regulated lncRNA-targets (d) in comparison of SD2_vs_NB. Pathways with P value < 0.05 were considered as significantly enriched

In addition, to further reveal the roles of DE lncRNAs, we performed Gene Ontology (GO) enrichment analysis to obtain the top 20 significantly enriched GO terms of biological processes, cellular components and molecular functions. With regard to the targets of up-regulated lncRNAs under SD, the most significantly enriched GO terms for biological processes were “pentose-phosphate shunt”, “isopentenyl diphosphate biosynthetic process, methylerythritol 4-phosphate pathway”, photosynthesis, light harvesting”, “nucleosome assembly”, “glucose catabolic process”, “rRNA processing”, “photosynthesis”, “cysteine biosynthetic process” and “tetrahydrofolate interconversion” (Fig. 5a). The GO terms including “chloroplast stroma”, “photosystem II oxygen evolving complex”, “nucleosome”, “thylakoid”, “chloroplast envelope”, “apoplast”, “chloroplast thylakoid”, “membrane” and “chloroplast” were enriched in cellular component category (Fig. 5a). The GO terms of “rRNA binding”, “enzyme inhibitor activity” and “cysteine synthase activity” were highly represented in molecular functions (Fig. 5a). The targets of the down-regulated lncRNAs under SD were abundant in 9, 6 and 5 GO terms of biological process, cellular component, and molecular function categories, respectively. As far as biological processes were concerned, “isopentenyl diphosphate biosynthetic process, methylerythritol 4-phosphate pathway”, “pentose-phosphate shunt”, “nucleosome assembly”, “reductive pentose-phosphate cycle”, “cysteine biosynthetic process”, “chloroplast relocation”, “thylakoid membrane organization”, “maltose metabolic process”, and “protein refolding” were the most important GO terms (Fig. 5b). The 6 GO terms such as “chloroplast stroma”, “chloroplast envelope”, “nucleosome”, “apoplast”, “chloroplast” and “thylakoid” were highlighted in cellular components (Fig. 5b). The GO terms like “rRNA binding”, “xyloglucan: xyloglucosyl transferase activity”, “protein heterodimerization activity”, “transferase activity, transferring glycosyl groups” and “L-ascorbate peroxidase activity” were predominant in molecular function category (Fig. 5b).

Fig. 5
figure5

Functional annotation of the differentially expressed lncRNAs by Gene Ontology (GO) classification. GO classifications of the up-regulated lncRNA-targets (a) and down-regulated lncRNA-targets (b) in comparison of SD1_vs_SD2. GO classifications of the up-regulated lncRNA-targets (c) and down-regulated lncRNA-targets (d) in comparison of SD2_vs_NB. The top 20 significant GO terms in biological process, cellular component and molecular function categories are selected based on the cutoff of P value < 0.05. Number right-side each horizonal column indicates the number of genes classified in the GO term

GO enrichment analysis of the targets of up-regulated lncRNAs under NB showed that “isopentenyl diphosphate biosynthetic process, methylerythritol 4-phosphate pathway”, “pentose-phosphate shunt”, “thylakoid membrane organization”, “nucleosome assembly”, “maltose metabolic process”, “chloroplast relocation”, “cysteine biosynthetic process”, “starch biosynthetic process”, “rRNA processing” and “protein refolding” were highly represented in biological processes (Fig. 5c). GO terms “chloroplast stroma”, “chloroplast envelope”, “thylakoid”, “nucleosome”, “chloroplast thylakoid membrane”, “chloroplast”, “apoplast” and “chloroplast thylakoid” showed high representativeness of cellular components (Fig. 5c). With regard to molecular function category, targets were involved in “rRNA binding” and “poly (U) RNA binding” (Fig. 5c). The targets of down-regulated lncRNAs under NB were significantly enriched in 9, 7 and 4 GO terms of biological process, cellular component, and molecular function categories, respectively. The biological process GO terms contained “pentose-phosphate shunt”, “isopentenyl diphosphate biosynthetic process, methylerythritol 4-phosphate pathway”, “nucleosome assembly”, “photosynthesis”, “photosynthesis, light harvesting”, “glucose catabolic process”, “cysteine biosynthetic process”, “cellular glucan metabolic process” and “chlorophyll biosynthetic process” (Fig. 5d). With regard to cellular components, “chloroplast stroma”, “thylakoid”, “nucleosome”, “chloroplast envelope”, “chloroplast thylakoid membrane”, “apoplast” and “photosystem II oxygen evolving complex” were highly represented (Fig. 5d). As far as molecular functions were concerned, “rRNA binding”, “cysteine synthase activity”, “xyloglucan: xyloglucosyl transferase activity” and “magnesium chelatase activity” were significantly enriched (Fig. 5d).

Prediction of the putative positive and negative flowering regulatory lncRNAs

According to the flowering time illustrated above, SD and NB exerted positive and negative effects on quinoa flowering, respectively. To ascertain which lncRNAs were probably involved in quinoa flowering, we analyzed the transcriptional profiles of the 91 DE lncRNAs. Through venn diagram analysis, we found that most of the DE lncRNAs (67 out of 91) were solely affected by SD or NB, and no DE lncRNAs were constantly induced or repressed by both SD and NB (Fig. 6a). Interestingly, we noticed that the expression patterns of 24 DE lncRNAs displayed opposite trends under SD compared with NB (Fig. 6a). Based on the expression profiles, we inferred that the 17 SD-induced but NB-repressed lncRNAs were the putative positive flowering regulator (Fig. 6b). On the other hand, the 7 lncRNAs that repressed by SD but activated by NB were predicted to be the putative flowering repressors (Fig. 6c).

Fig. 6
figure6

Identification and expression analysis of the putative positive and negative flowering lncRNAs. a Venn diagram showing the number of specific and common differentially expressed lncRNAs in comparisons of SD1_vs_SD2 and SD2_vs_NB. The 17 DE lncRNAs showing up-regulation in SD1_vs_SD2 but down-regulation in SD2_vs_NB were predicted to be putative positive flowering lncRNAs. The 7 DE lncRNAs showing down-regulation in SD1_vs_SD2 but up-regulation in SD2_vs_NB were predicted to be putative negative flowering lncRNAs. Heatmap showing the expression levels (log10FPKM) of putative positive (b) and negative (c) flowering lncRNAs in different samples

Further, to understand the functions of the lncRNAs, we investigated the regulatory networks of the putative flowering lncRNAs. With regard to the 17 positive flowering lncRNAs, 143 cis targets with 147 cis-acting lncRNA-mRNA pairs and 4034 trans targets with 16,844 trans-acting lncRNA-mRNA pairs were selected, respectively (Fig. 7) (Additional file 8: Table S7). Cis-acting network analysis indicated that the largest node MSTRG.39390.2 co-located with 19 targets and MSTRG.59128.3 only co-located with 1 target (Fig. 7a-b). Notably, a few important plant growth and development related transcription factors (TFs) were identified among the cis targets. For example, 3 SCARECROW-like genes of GRAS family (AUR62001765, AUR62001766, AUR62001768) co-located with MSTRG.39390.2; 1 WOX (WUSCHEL-related homeobox) (AUR62035466) co-located with MSTRG.47393.1. According to the trans-acting network of the positive flowering lncRNAs, the largest node MSTRG.39390.2 and smallest node MSTRG.10578.1 were predicted to regulate 1495 and 343 trans targets, respectively (Fig. 7c-d). Moreover, we found that 1 mRNA could be the common target of 1 to 14 lncRNAs and 423 mRNA genes were correlated with more than 10 lncRNAs (Fig. 7d). The expression heatmap demonstrated that most of these trans targets harbored the same expression trend with the lncRNAs (Fig. 7c). Functional annotation by KEGG enrichment analysis suggested that no pathways were significantly enriched for cis targets, whereas the trans target genes were highly abundant in pathways “Ribosome biogenesis in eukaryotes”, “Photosynthesis-antenna proteins” and “Photosynthesis” (Table 3).

Fig. 7
figure7

Expression and regulatory network analysis of the cis and trans targets of positive flowering lncRNAs. a Heatmap showing the expression levels (log10FPKM) of the cis targets in different samples. b Relationship between the positive flowering lncRNAs and their co-located cis targets. c Heatmap showing the expression levels (log10FPKM) of the trans targets in different samples. d Co-expression network analysis of the positive flowering lncRNAs and their trans targets. The orange triangles and light-blue circles stand for lncRNAs and mRNA targets, respectively. The node and font size reflect the degrees between lncRNAs and their targets

Table 3 KEGG enrichment analysis of the cis and trans targets of the putative positive and negative flowering lncRNAs

For the 7 putative flowering repressor lncRNAs, 50 cis targets with 57 cis-acting lncRNA-mRNA pairs and 2615 trans targets with 3985 trans-acting lncRNA-mRNA pairs may be correlated with these lncRNAs (Fig. 8) (Additional file 9: Table S8). From the cis-acting network we found that the largest node MSTRG.5158.2 co-located with 18 targets, while the smallest nodes MSTRG.53750.2 and MSTRG.3121.1 only co-located with 1 target (Fig. 8a-b). With regard to the trans co-expression network, the largest node MSTRG.29609.1 and the smallest node MSTRG.5158.2 may interact with 1164 and 44 targets, respectively, and 1 mRNA could be correlated with 1 to 4 lncRNAs (Fig. 8c-d). The expression heatmap showed that about 50% of these trans targets harbored similar expression trend with the negative flowering lncRNAs (Fig. 8c). KEGG enrichment analysis indicated the cis targets were abundant in pathway “Glutathione metabolism”, and the trans targets were significantly enriched in pathways “Photosynthesis-antenna proteins”, “Ribosome biogenesis in eukaryotes”, “Porphyrin and chlorophyll metabolism” and “Glyoxylate and dicarboxylate metabolism” (Table 3).

Fig. 8
figure8

Expression and regulatory network analysis of the cis and trans targets of negative flowering lncRNAs. a Heatmap showing the expression levels (log10FPKM) of the cis targets in different samples. b Relationship between the negative flowering lncRNAs and their co-located cis targets. c Heatmap showing the expression levels (log10FPKM) of the trans targets in different samples. d Co-expression network analysis of the negative flowering lncRNAs and their trans targets. The orange triangles and light-blue circles stand for lncRNAs and mRNA targets, respectively. The node and font size reflect the degrees between lncRNAs and their targets

Co-expressed network associated with quinoa flowering

Furthermore, by comparing the seed nodes in the positive and negative lncRNAs flowering networks, we noticed that 1653 mRNAs were simultaneously regulated by both types of regulatory lncRNAs, accounting for 39.66% (1653/4168) and 62.12% (1653/2661) of the positive and negative network targets, respectively (Fig. 9a, b). The expression heatmap indicated that most of these targets harbor higher expression levels under SD than under NB (Fig. 9c). Co-expression network analysis demonstrated that the largest node MSTRG.69967.2 was linked with 1176 targets while MSTRG.59128.3 was only co-expressed with 1 target, and the top-5-degree target nodes (Chenopodium_quinoa_newGene_17423, AUR62040917, AUR62018533, AUR62015897, AUR62000858) were regulated by as many as 17 lncRNAs (Fig. 9b). Intriguingly, among the 1653 targets, we found that a few vital plant growth and development related modulators co-expressed with 15 positive and 5 negative putative flowering lncRNAs (Fig. 10a-b). For example, the FLOWERING LOCUS T-like (FT-like) homolog (AUR62006619), which belongs to the PHOSPHATIDYLETHANOLAMINE-BINDING PROTEIN (PEBP) family and functions as a crucial floral integrator triggering downstream floral organogenesis in plants [10], was predicted to be co-expressed with 13 positive (MSTRG.53484.1, MSTRG.30736.1, MSTRG.39390.2, MSTRG.27734.1, MSTRG.14725.1, MSTRG.47393.1, MSTRG.3107.2, MSTRG.29153.1, MSTRG.24870.1, MSTRG.10889.3, MSTRG.8116.1, MSTRG.69967.1, MSTRG.69967.2) and 1 negative (MSTRG.29609.1) flowering lncRNAs (Fig. 10a-b). Another PEBP member TWIN SISTER of FT-like (TSF-like) homolog (AUR62026433), which shares similar functions with FT, may be regulated by 5 positive (MSTRG.30736.1, MSTRG.39390.2, MSTRG.47393.1, MSTRG.28884.1, MSTRG.69967.2) and 2 negative (MSTRG.3121.1, MSTRG.29609.1) flowering lncRNAs (Fig. 10a-b). AUR62039984, encoding a B-BOX Zinc Finger protein CONSTANS-like (CO-like) which acts upstream of FT in photoperiod-dependent manner [44], was tightly correlated with the flowering activator MSTRG.8116.1 and the flowering repressor MSTRG.32307.3. AUR62009205, an EARLY FLOWERING 3-like homolog which plays important roles in maintaining circadian rhythms and controlling flowering time [45, 46], was the trans target of 2 positive (MSTRG.32307.3, MSTRG.3121.1) and 2 negative (MSTRG.39390.2, MSTRG.69967.1) flowering lncRNAs (Fig. 10a-b). The MYB-type LATE ELONGATED HYPOCOTYL (LHY) is a core component of circadian clock controlling many clock outputs such as flowering and photomorphogenesis [47]. In this study, we found that the quinoa LHY homolog (AUR62004570) was co-expressed with 6 lncRNAs (5 positive lncRNAs: MSTRG.39390.2, MSTRG.27734.1, MSTRG.29153.1, MSTRG.8116.1, MSTRG.69967.1; 1 negative lncRNA: MSTRG.3121.1) (Fig. 10a-b). The bZIP transcription factor ELONGATED HYPOCOTYL 5 (HY5) mediates blue light signaling to modulate circadian clock rhythms and plays multifaceted roles in plant growth and development [48]. The quinoa HY5 homolog AUR62030640 was correlated with 4 positive (MSTRG.39390.2, MSTRG.29153.1, MSTRG.8116.1, MSTRG.69967.2) and 2 negative lncRNAs (MSTRG.32307.3, MSTRG.29609.1) (Fig. 10a-b). Previous studies suggested that the RELATED TO ABI3/VP1 (RAV) family members delay heading time by repressing Heading date 3a (Hd3a) and GIBBERELLIN 3-OXIDASE1/2 in rice (Oryza sativa) [49]. In the present study, we found that the quinoa RAV1 homolog (AUR62035279) was regulated by 3 lncRNAs (2 positive lncRNAs, MSTRG.39390.2, MSTRG.69967.1; 1 negative lncRNA: MSTRG.3121.1) (Fig. 10a-b). Numerous evidences have indicated that photoreceptor PHYTOCHROME A (PHYA) mainly mediates light responses to various ratios of red/far-red light. In rice, the phyA mutation in phyB or phyC background leads to dramatic early heading [14]. The quinoa PHYA homolog (AUR62003557) was identified as the target of 8 lncRNAs (5 positive lncRNAs, MSTRG.30736.1, MSTRG.29153.1, MSTRG.10889.3, MSTRG.10578.1, MSTRG.69967.2; 3 negative lncRNAs: MSTRG.11610.1, MSTRG.29609.4, MSTRG.29609.1) (Fig. 10a-b). In addition, we found another photoreceptor, the blue light photoreceptor CRYPTOCHROME1 homolog (AUR62018922) which is involved in entrainment of circadian clock and photoperiodic flowering [50], was linked with 7 lncRNAs (4 positive lncRNAs, MSTRG.30736.1, MSTRG.47393.1, MSTRG.28884.1, MSTRG.69967.2; 3 negative lncRNAs: MSTRG.11610.1, MSTRG.29609.4, MSTRG.29609.1) (Fig. 10a-b). Consequently, the tight association with these photoreceptors, circadian clock and florigen genes indicates these lncRNAs may play important roles in photoperiodic flowering of quinoa.

Fig. 9
figure9

Identification of the common targets of positive and negative flowering lncRNAs. a Venn diagram showing the number of common and specific targets of positive and negative flowering lncRNAs. b Co-expression network of 17 positive, 7 negative flowering lncRNAs and their common targets. The orange triangles, green V and light-blue circles represent positive, negative flowering lncRNAs and their common targets, respectively. The node and font size reflect the degrees between lncRNAs and their targets. c Heatmap showing the expression levels (log10FPKM) of the 1653 common targets in different samples

Fig. 10
figure10

Construction of the core network involved in quinoa flowering. a Co-expression network of the positive, negative flowering lncRNAs and vital flowering homologs. The orange triangles, green rhombuses and purple circles represent positive, negative flowering lncRNAs and the vital flowering homologs, respectively. b Heatmap showing the expression levels (log10FPKM) of the vital flowering homologs in different samples

Discussion

Flowering is a vital process closely associated with the final yield of crops. Previous studies have revealed the roles of many protein-coding genes in flowering time regulation. However, knowledge about the functions of lncRNAs in plant flowering and yield regulation is still limited. Until recently, a few studies have deciphered the pivotal roles of lncRNAs in flowering and yield regulation. In Arabidopsis, lncRNAs COLDAIR and COOLAIR were proven to regulate vernalization and flowering time via epigenic-repression of FLC [23, 24]. In rice, EARLY FLOWERING-COMPLETELY DOMINANT (EF-CD) encodes an antisense lncRNA overlapping with the floral activator OsSOC1 to shorten maturation period without concomitant yield penalty [51]. Mining of specific lncRNAs may provide novel strategies for crop yield improvement.

Quinoa is an annual-flowering crop originated from the low latitude Andeans region. The flowering time and yield of some quinoa cultivars tend to vary when introduced into different countries [52]. However, the flowering regulators, especially the lncRNAs mediated flowering networks, remain to be revealed. SD combined with NB approach is an effective way for exploring the floral regulators. In this study, we found NB repressed the quinoa floral buds development in contrast to the flowering-promoting effects of SD (Fig. 1a-d). These results were consistent with the observations in other SDPs such as rice [11] and soybean [12], indicating the conserved flowering-repressing effects of NB on SDPs. To uncover the specific lncRNAs involved in NB responses, we compared the transcriptome differences between plants grown under SD and NB. By high-throughput sequencing, we identified 4914 lncRNAs from quinoa leaf tissues, less than the number of lncRNAs identified from tomato (Solanum lycopersicum) [53] but more than that identified from Cucumis melo [31], Cassava (Manihot esculenta) [54] and peanut (Arachis hypogaea) [55]. Meanwhile, only 605 (12.3%) shared similar sequences with the known lncRNAs in the NONCODE database. We presumed that the lncRNAs number and sequence conservation may be highly associated with tissue specificity and plant species. In this study, we identified 91 DE lncRNAs, most of which (73.6%) were affected solely either by SD or NB. Considering the contrary effects of SD and NB, we neglected this part of DE lncRNAs, and speculated that 17 and 7 lncRNAs, affected by both SD and NB in antagonizing manner, were the putative positive and negative flowering regulators, respectively (Fig. 6a-c).

To investigate the roles of the lncRNAs, we performed KEGG functional enrichment analysis of the cis and trans targets of the putative flowering lncRNAs. Pathway “Glutathione metabolism” was specific for the cis targets of negative flowering lncRNAs. “Photosynthesis” was the specific pathway for the trans targets of positive flowering lncRNAs, while “Porphyrin and chlorophyll metabolism” and “Glyoxylate and dicarboxylate metabolism” were specific for the trans targets of negative flowering lncRNAs (Table 3). These specific pathways may provide valuable cues for investigating the roles of these flowering lncRNAs in future.

Comparison of the trans targets of positive and negative flowering lncRNAs revealed that a considerable number of targets (1653 genes) co-expressed with both types of lncRNAs (Fig. 9a), indicative of competing flowering regulatory networks in quinoa. As supporting evidences, annotation of the common targets revealed a few floral homologs which may function as vital mediators of the antagonizing lncRNAs. Flowering time is largely influenced by endogenous circadian clock and external light conditions which are output to the final integrators FT in leaves [4]. As discovered in a lot of plant species, there are multiple FT homologs with highly conserved sequences in their genomes. However, different FT homologs may experience functional diversification because of the varied expression patterns [56]. For example, in rice there are 11 FT homologous genes, among which both Hd3a and RFT1 are floral activators [11, 57]. However, only Hd3a showed obvious expression suppression by NB [11], which is the principal cause of the NB effect in rice. Ten FT homologs are isolated in soybean genome, yet, according to the expression pattern they are divided into SD-, LD-specific and photoperiod-independent groups [58]. Functional characterization suggests that GmFT1a and GmFT2a/5a have opposite effects in flowering of soybean [58]. In Chenopodium rubrum (C. rubrum), a close relative of quinoa, CrFTL1 is SD-inducible and inhibited by NB, whereas CrFTL2 is constitutively expressed [59]. In allotetraploid quinoa, 11 putative FT homologs have been identified [60]. However, which are NB-responsive is not determined. In this study, we found that only AUR62006619 (FT-like) and AUR62026433 (TSF-like) were significantly affected by NB (Fig. 10). AUR62006619 and AUR62026433 were down-regulated by 8-fold and 2-fold, respectively, under NB than under SD. Thus, we speculated that these two FT homologs are probably the essential mediators for NB responses in quinoa. Co-expression network analysis showed that the quinoa FT-like homolog (AUR62006619) and TSF-like homolog (AUR62026433) were regulated by 14 and 7 lncRNAs, respectively. Noteworthy, both of these two quinoa FT homologs were commonly regulated by 5 lncRNAs (4 positive lncRNAs: MSTRG.30736.1, MSTRG.39390.2, MSTRG.47393.1, MSTRG.69967.2; 1 negative lncRNAs: MSTRG.29609.1) (Fig. 10), suggesting these two homologs may function through partially overlapped signal cascade.

The B-BOX type CO, a direct activator of FT, functions as the output of circadian clock [61]. In rice, the transcriptional level of CO homolog-HD1 is not obviously changed by NB [11]. However, different scenarios are observed in other plants. For instance, in Chrysanthemum (Chrysanthemum morifolium), the expression of CmCOL1 is slightly reduced by NB [62]. In C. rubrum, two CrCOLs genes were down-regulated by dark-light transition [59]. Similarly, we found that 1 out 6 CO-like homologs (AUR62039984) in quinoa was down-regulated by NB. These results indicate the specific roles of CO homologs of different plants in NB responses. Coincidently, three circadian clock related genes, LHY-like (AUR62004570), ELF3-like (AUR62009205) and HY5-like (AUR62030640), were also repressed by NB. In addition, LHY-like, HY5-like and CO-like were commonly regulated by one lncRNA (MSTRG.8116.1) (Fig. 10). These results imply the possible positive roles of these circadian clock genes on quinoa CO-like. Photoreceptors are essential for perceiving light signals to regulate circadian clock and flowering time. Phytochromes and Cryptochromes are responsible for red/far-red and blue light perceiving, respectively. In quinoa, there are 6 Phytochromes (2 PHYA, 2 PHYB and 2 PHYC) and 4 Cryptochromes (2 CRY1 and 2 CRY2) [60]. In this study, we found 1 PHYA and 1 CRY1 were specifically up-regulated by NB and were commonly co-expressed with as many as 5 lncRNAs (MSTRG.69967.2, MSTRG.30736.1, MSTRG.29609.1, MSTRG.29609.4, MSTRG.11610.1 and MSTRG.30736.1) (Fig. 10). It is of note that MSTRG.69967.2, MSTRG.30736.1 and MSTRG.29609.1 also participated in regulating the two FT homologs (Fig. 10). Thus, these data suggested PHYA and CRY1 may regulate flowering time by controlling the two FT homologs. Together, these putative flowering lncRNAs may mediate NB responses by regulating the key flowering genes. We believe further functional investigation on these lncRNAs will shed light on their roles in quinoa flowering.

Conclusion

In summary, in this study we performed SS RNA-seq to explore genome-wide lncRNAs in quinoa. By analyzing the transcriptome data under SD and NB, we predict the specific lncRNAs that may be required for photoperiodic flowering. To our knowledge, this is the first systematic investigation on the lncRNAs involved in NB-responses of quinoa. In total, we identified 17 putative positive- and 7 putative negative-flowering lncRNAs. Co-expression network analysis indicated that 15 positive-flowering lncRNAs may compete with 5 negative-flowering lncRNAs to control a few pivotal flowering homologous genes, which might be the basis of NB responses in quinoa. Collectively, these results will deepen our understanding of the roles of lncRNAs in flowering, and provide valuable genetic resources for yield improvement of quinoa.

Methods

Plant materials and growth conditions

Quinoa (Chenopodium quinoa) cultivar “HL1” and “HL2” were used for flowering time analysis. The two cultivars have nearly the same floral transition time, however, “HL2” needs longer time (about 115 d) compared with “HL1” (about 60 d) for post-flowering development. Seeds were harvested from the Experimental Station of Key Laboratory of Coarse Cereal Processing, Ministry of Agriculture and Rural Affairs, in 2019 at Liangshan, Sichuan province. Seeds were sown in a pot with 20 cm in diameter and 18 cm in depth. Quinoa seedlings were maintained for 14 d in a greenhouse at Chengdu University under natural day conditions. Then, the seedlings pots were transferred into growth chambers (Jiang Nan, Ning Bo, China) subjecting to SD and NB treatment. The floral buds developmental situation of both cultivers under different conditions was observed every day.

For RNA-seq, the seeds of “HL1” were sown in pots with the same size as described above. After sowing, the pots were placed in growth chambers for SD treatment with daily cycles of 8 h light at 25 °C and 16 h dark at 23 °C. The relative humidity was 70%, and light intensity was 750 μmol.m− 2.s− 1. On 7 DAS, 15 robust shoots were retained in each pot by removing the spared ones. On 14 DAS (SD1) and 26 DAS (SD2), the top two fully expanded leaves of each seedling were collected at 2 h after dawn, respectively. Then, the growth conditions were reset to night-break conditions (NB) with daily cycles of 8 h light at 25 °C and 7 h dark, 1 h light exposure and 8 h dark at 23 °C. After 2 d NB treatment, namely on 28 DAS, the leaf samples were collected as above. Sampling was performed with three replicates and each replicate contains 5 individual plants. The harvested fresh leaf samples were immediately frozen by liquid nitrogen and stored at − 80 °C before RNA extraction.

Construction of strand-specific RNA sequencing libraries

Total RNA was extracted from each leaves sample using Trizol reagent (Invitrogen, CA, USA) according to the manufacturer’s instructions. RNA degradation was checked by 1.5% agarose gel electrophoresis. RNA concentration and purity of each sample were measured to be qualified for downstream experiments. Only the samples with RIN (RNA Integrity Number) value larger than 7.5 were used for long-noncoding RNA sequencing libraries construction. An amount of 1.5 μg total RNA per sample was input for rRNA removal reaction by using the Ribo-Zero rRNA Removal Kit (Epicentre, Madison, WI, USA). Then, strand-specific RNA sequencing libraries were generated using NEBNextR Ultra™ Directional RNA Library Prep Kit for IlluminaR (NEB, USA) according to manufacturer’s protocol.

Data analysis

RNA sequencing was performed on an Illumina platform (Hiseq 4000) to generate 150 bp paired-end reads. Raw reads from 9 libraries were filtered by removing the unqualified reads (adapter-, ploy-N-containing and low-quality reads) to obtain clean reads. By applying HISAT2 software [38], clean reads were mapped to the quinoa reference genome sequences (https://www.cbrc.kaust.edu.sa/chenopodiumdb) [37]. Based on the mapping results, StringTie [39] was used to re-construct the transcripts. Gffcompare program [40] was used to annotate the assembled transcripts, and the unannotated transcripts were further used for long-noncoding RNA prediction. Transcripts with lengths > 200 bp and exon number ≥ 2 were retained as lncRNA candidates whose protein-coding potential were further evaluated by using CPC (Coding Potential Calculator) [63], CNCI (Coding-Non-Coding Index) [64], CPAT (Coding Potential Assessment Tool) [65] and Pfam programs [66] to obtain the final lncRNAs. LncRNAs were categorized into different types of lncRNAs including lincRNA, intronic lncRNA, anti-sense lncRNA and sense lncRNA by using Cuffcompare program.

Differential expression analysis

The expression levels of both mRNA and lncRNA genes were calculated and normalized to FPKM (fragments per kilo-base of exon per million fragments mapped) [42] value using StringTie [39]. Differentially expressed genes (DEGs) between different groups were determined by using DESeq2 R package [43] based on the negative binomial distribution model. The resulting P values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate (FDR). Thresholds of an FDR < 0.05 and |log2 (fold change)| ≥ 1 were applied to call DEGs.

Target gene prediction

For cis targets identification, the protein-coding genes located at 100 kb upstream and downstream of lncRNA were selected and functionally annotated. To identify the trans targets, the correlation between mRNAs and lncRNAs was evaluated using the Pearson’s correlation coefficient (PCC) based on their expression profiles at different samples. mRNA-lncRNA pairs with |PCC value| > 0.9 and P value < 0.01 were designated as co-expressed gene pairs in which the mRNA genes were identified as the trans targets of lncRNAs.

Functional enrichment analysis of the targets of differentially expressed lncRNAs

To understand the biological function of differentially expressed lncRNAs (DE lncRNAs), Kyoto Encyclopedia of Genes and Genomes (KEGG) [67] and Gene Ontology enrichment analysis [68] of the corresponding target genes was conducted. KEGG enrichment analysis was performed by using KOBAS software and the pathways with corrected P value less than 0.05 were taken as significantly enriched by the target genes. GO enrichment analysis of the target genes was implemented by the topGO R packages. The most significantly enriched GO terms were selected with the cutoff of corrected P value < 0.05.

Co-expression network analysis and visualization

Based on the expression data, the PCC between lncRNAs and mRNAs was obtained. Co-expressed lncRNA-mRNA pairs were visualized using Cytoscape (Version 3.7.2) [69]. By using the Cytoscape tools “Network Analyzer-Network analysis”, the network was analyzed and the degrees of nodes were obtained.

Availability of data and materials

The transcriptome raw data produced in this study is deposited at NCBI SRA database (http://trace.ncbi.nlm.nih.gov/Traces/sra) under accession PRJNA673052.

Abbreviations

SD:

Short-day

NB:

Night-break

lncRNAs:

Long noncoding RNAs

FPKM:

Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced

SDP:

Short-day plant

LDP:

Long-day plant

DAS:

Days after sowing

GO:

Gene ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

PCC:

Pearson’s correlation coefficient

CPC:

Coding Potential Calculator

CNCI:

Coding-Non-Coding Index

CPAT:

Coding Potential Assessment Tool

References

  1. 1.

    Fornara F, de Montaigu A, Coupland G. SnapShot: Control of flowering in Arabidopsis. Cell. 2010;141(3):550, 550.e551–2.

    Article  Google Scholar 

  2. 2.

    Wigge PA. FT, a mobile developmental signal in plants. Curr Biol. 2011;21(9):R374–8. https://doi.org/10.1016/j.cub.2011.03.038.

    CAS  Article  PubMed  Google Scholar 

  3. 3.

    Andres F, Coupland G. The genetic basis of flowering responses to seasonal cues. Nat Rev Genet. 2012;13(9):627–39. https://doi.org/10.1038/nrg3291.

    CAS  Article  PubMed  Google Scholar 

  4. 4.

    Song YH, Ito S, Imaizumi T. Flowering time regulation: photoperiod- and temperature-sensing in leaves. Trends Plant Sci. 2013;18(10):575–83. https://doi.org/10.1016/j.tplants.2013.05.003.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  5. 5.

    Taoka K, Ohki I, Tsuji H, Furuita K, Hayashi K, Yanase T, et al. 14-3-3 proteins act as intracellular receptors for rice Hd3a florigen. Nature. 2011;476(7360):332–5. https://doi.org/10.1038/nature10272.

    CAS  Article  PubMed  Google Scholar 

  6. 6.

    Putterill J, Varkonyi-Gasic E. FT and florigen long-distance flowering control in plants. Curr Opin Plant Biol. 2016;33:77–82. https://doi.org/10.1016/j.pbi.2016.06.008.

    CAS  Article  PubMed  Google Scholar 

  7. 7.

    Itoh H, Nonoue Y, Yano M, Izawa T. A pair of floral regulators sets critical day length for Hd3a florigen expression in rice. Nat Genet. 2010;42(7):635–8. https://doi.org/10.1038/ng.606.

    CAS  Article  PubMed  Google Scholar 

  8. 8.

    Kong F, Liu B, Xia Z, Sato S, Kim BM, Watanabe S, et al. Two coordinately regulated homologs of FLOWERING LOCUS T are involved in the control of photoperiodic flowering in soybean. Plant Physiol. 2010;154(3):1220–31. https://doi.org/10.1104/pp.110.160796.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  9. 9.

    Shimada S, Ogawa T, Kitagawa S, Suzuki T, Ikari C, Shitsukawa N, et al. A genetic network of flowering-time genes in wheat leaves, in which an APETALA1/FRUITFULL-like gene, VRN1, is upstream of FLOWERING LOCUS T. Plant J. 2009;58(4):668–81. https://doi.org/10.1111/j.1365-313X.2009.03806.x.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  10. 10.

    Corbesier L, Vincent C, Jang S, Fornara F, Fan Q, Searle I, et al. FT protein movement contributes to long-distance signaling in floral induction of Arabidopsis. Science. 2007;316(5827):1030–3. https://doi.org/10.1126/science.1141752.

    CAS  Article  PubMed  Google Scholar 

  11. 11.

    Ishikawa R, Tamaki S, Yokoi S, Inagaki N, Shinomura T, Takano M, et al. Suppression of the floral activator Hd3a is the principal cause of the night break effect in rice. Plant Cell. 2005;17(12):3326–36. https://doi.org/10.1105/tpc.105.037028.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  12. 12.

    Xu M, Yamagishi N, Zhao C, Takeshima R, Kasai M, Watanabe S, et al. The soybean-specific maturity gene E1 family of floral repressors controls night-break responses through Down-regulation of FLOWERING LOCUS T Orthologs. Plant Physiol. 2015;168(4):1735–46. https://doi.org/10.1104/pp.15.00763.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  13. 13.

    Pearce S, Shaw LM, Lin H, Cotter JD, Li C, Dubcovsky J. Night-break experiments shed light on the Photoperiod1-mediated flowering. Plant Physiol. 2017;174(2):1139–50. https://doi.org/10.1104/pp.17.00361.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  14. 14.

    Ishikawa R, Shinomura T, Takano M, Shimamoto K. Phytochrome dependent quantitative control of Hd3a transcription is the basis of the night break effect in rice flowering. Genes Genet Syst. 2009;84(2):179–84. https://doi.org/10.1266/ggs.84.179.

    CAS  Article  PubMed  Google Scholar 

  15. 15.

    Ramos-Sánchez JM, Triozzi PM, Alique D, Geng F, Gao M, Jaeger KE, et al. LHY2 Integrates Night-Length Information to Determine Timing of Poplar Photoperiodic Growth. Curr Biol. 2019;29(14):2402–2406.e2404.

    Article  Google Scholar 

  16. 16.

    Chekanova JA. Long non-coding RNAs and their functions in plants. Curr Opin Plant Biol. 2015;27:207–16. https://doi.org/10.1016/j.pbi.2015.08.003.

    CAS  Article  PubMed  Google Scholar 

  17. 17.

    Liu X, Hao L, Li D, Zhu L, Hu S. Long non-coding RNAs and their biological roles in plants. Genomics Proteomics Bioinformatics. 2015;13(3):137–47. https://doi.org/10.1016/j.gpb.2015.02.003.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  18. 18.

    Hu S, Wang X, Shan G. Insertion of an Alu element in a lncRNA leads to primate-specific modulation of alternative splicing. Nat Struct Mol Biol. 2016;23(11):1011–9. https://doi.org/10.1038/nsmb.3302.

    CAS  Article  PubMed  Google Scholar 

  19. 19.

    Wei S, Chen H, Dzakah EE, Yu B, Wang X, Fu T, et al. Systematic evaluation of C. elegans lincRNAs with CRISPR knockout mutants. Genome Biol. 2019;20(1):7.

    Article  Google Scholar 

  20. 20.

    Rashid F, Shah A, Shan G. Long non-coding RNAs in the cytoplasm. Genomics Proteomics Bioinformatics. 2016;14(2):73–80. https://doi.org/10.1016/j.gpb.2016.03.005.

    Article  PubMed  PubMed Central  Google Scholar 

  21. 21.

    Yu B, Shan G. Functions of long noncoding RNAs in the nucleus. Nucleus. 2016;7(2):155–66. https://doi.org/10.1080/19491034.2016.1179408.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  22. 22.

    Kim DH, Xi Y, Sung S. Modular function of long noncoding RNA, COLDAIR, in the vernalization response. PLoS Genet. 2017;13(7):e1006939. https://doi.org/10.1371/journal.pgen.1006939.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  23. 23.

    Heo JB, Sung S. Vernalization-mediated epigenetic silencing by a long intronic noncoding RNA. Science. 2011;331(6013):76–9. https://doi.org/10.1126/science.1197349.

    CAS  Article  PubMed  Google Scholar 

  24. 24.

    Sun Q, Csorba T, Skourti-Stathaki K, Proudfoot NJ, Dean C. R-loop stabilization represses antisense transcription at the Arabidopsis FLC locus. Science. 2013;340(6132):619–21. https://doi.org/10.1126/science.1234848.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  25. 25.

    Ding J, Lu Q, Ouyang Y, Mao H, Zhang P, Yao J, et al. A long noncoding RNA regulates photoperiod-sensitive male sterility, an essential component of hybrid rice. Proc Natl Acad Sci U S A. 2012;109(7):2654–9. https://doi.org/10.1073/pnas.1121374109.

    Article  PubMed  PubMed Central  Google Scholar 

  26. 26.

    Franco-Zorrilla JM, Valli A, Todesco M, Mateos I, Puga MI, Rubio-Somoza I, et al. Target mimicry provides a new mechanism for regulation of microRNA activity. Nat Genet. 2007;39(8):1033–7. https://doi.org/10.1038/ng2079.

    CAS  Article  PubMed  Google Scholar 

  27. 27.

    Liu X, Li D, Zhang D, Yin D, Zhao Y, Ji C, et al. A novel antisense long noncoding RNA, TWISTED LEAF, maintains LEAF blade flattening by regulating its associated sense R2R3-MYB gene in rice. New Phytol. 2018;218(2):774–88. https://doi.org/10.1111/nph.15023.

    CAS  Article  PubMed  Google Scholar 

  28. 28.

    Campalans A, Kondorosi A, Crespi M. Enod40, a short open reading frame-containing mRNA, induces cytoplasmic localization of a nuclear RNA binding protein in Medicago truncatula. Plant Cell. 2004;16(4):1047–59. https://doi.org/10.1105/tpc.019406.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  29. 29.

    Ariel F, Jegu T, Latrasse D, Romero-Barrios N, Christ A, Benhamed M, et al. Noncoding transcription by alternative RNA polymerases dynamically regulates an auxin-driven chromatin loop. Mol Cell. 2014;55(3):383–96. https://doi.org/10.1016/j.molcel.2014.06.011.

    CAS  Article  PubMed  Google Scholar 

  30. 30.

    Huang L, Dong H, Zhou D, Li M, Liu Y, Zhang F, et al. Systematic identification of long non-coding RNAs during pollen development and fertilization in Brassica rapa. Plant J. 2018;96(1):203–22. https://doi.org/10.1111/tpj.14016.

    CAS  Article  PubMed  Google Scholar 

  31. 31.

    Tian Y, Bai S, Dang Z, Hao J, Zhang J, Hasi A. Genome-wide identification and characterization of long non-coding RNAs involved in fruit ripening and the climacteric in Cucumis melo. BMC Plant Biol. 2019;19(1):369. https://doi.org/10.1186/s12870-019-1942-4.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  32. 32.

    Fu L, Ding Z, Tan D, Han B, Sun X, Zhang J. Genome-wide discovery and functional prediction of salt-responsive lncRNAs in duckweed. BMC Genomics. 2020;21(1):212.

    CAS  Article  Google Scholar 

  33. 33.

    Tan X, Li S, Hu L, Zhang C. Genome-wide analysis of long non-coding RNAs (lncRNAs) in two contrasting rapeseed (Brassica napus L.) genotypes subjected to drought stress and re-watering. BMC Plant Biol. 2020;20(1):81.

    CAS  Article  Google Scholar 

  34. 34.

    Lopez-Marques RL, Norrevang AF, Ache P, Moog M, Visintainer D, Wendt T, et al. Prospects for the accelerated improvement of the resilient crop quinoa. J Exp Bot. 2020;71(18):5333–47. https://doi.org/10.1093/jxb/eraa285.

    Article  PubMed  PubMed Central  Google Scholar 

  35. 35.

    Storchova H, Hubackova H, Abeyawardana OAJ, Walterova J, Vondrakova Z, Eliasova K, et al. Chenopodium ficifolium flowers under long days without upregulation of FLOWERING LOCUS T (FT) homologs. Planta. 2019;250(6):2111–25. https://doi.org/10.1007/s00425-019-03285-1.

    CAS  Article  PubMed  Google Scholar 

  36. 36.

    Christiansen JL, Jacobsen SE, Jørgensen ST. Photoperiodic effect on flowering and seed development in quinoa (Chenopodium quinoaWilld.). Acta Agric Scand B Soil Plant Sci. 2010;60(6):539–44. https://doi.org/10.1080/09064710903295184.

    Article  Google Scholar 

  37. 37.

    Jarvis DE, Ho YS, Lightfoot DJ, Schmockel SM, Li B, Borm TJ, et al. The genome of Chenopodium quinoa. Nature. 2017;542(7641):307–12. https://doi.org/10.1038/nature21370.

    CAS  Article  PubMed  Google Scholar 

  38. 38.

    Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60. https://doi.org/10.1038/nmeth.3317.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  39. 39.

    Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5. https://doi.org/10.1038/nbt.3122.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  40. 40.

    Pertea G, Pertea M. GFF Utilities: GffRead and GffCompare. F1000Res. 2020;9:ISCB Comm J-304.

    Article  Google Scholar 

  41. 41.

    Liu C, Bai B, Skogerbo G, Cai L, Deng W, Zhang Y, et al. NONCODE: an integrated knowledge database of non-coding RNAs. Nucleic Acids Res. 2005;33(Database issue):D112–5. https://doi.org/10.1093/nar/gki041.

    CAS  Article  PubMed  Google Scholar 

  42. 42.

    Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5. https://doi.org/10.1038/nbt.1621.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  43. 43.

    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. https://doi.org/10.1186/s13059-014-0550-8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  44. 44.

    Luccioni L, Krzymuski M, Sanchez-Lamas M, Karayekov E, Cerdan PD, Casal JJ. CONSTANS delays Arabidopsis flowering under short days. Plant J. 2019;97(5):923–32. https://doi.org/10.1111/tpj.14171.

    CAS  Article  PubMed  Google Scholar 

  45. 45.

    Yang Y, Peng Q, Chen GX, Li XH, Wu CY. OsELF3 is involved in circadian clock regulation for promoting flowering under long-day conditions in rice. Mol Plant. 2013;6:202–15. https://doi.org/10.1093/mp/sss062.

  46. 46.

    Zhao J, Huang X, Ouyang X, Chen W, Du A, Zhu L, et al. OsELF3-1, an Ortholog of Arabidopsis EARLY FLOWERING 3, regulates Rice circadian rhythm and photoperiodic flowering. PLoS One. 2012;7(8):e43705. https://doi.org/10.1371/journal.pone.0043705.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  47. 47.

    Alabadi D, Yanovsky MJ, Mas P, Harmer SL, Kay SA. Critical role for CCA1 and LHY in maintaining circadian rhythmicity in Arabidopsis. Curr Biol. 2002;12(9):757–61. https://doi.org/10.1016/S0960-9822(02)00815-1.

    CAS  Article  PubMed  Google Scholar 

  48. 48.

    Hajdu A, Dobos O, Domijan M, Balint B, Nagy I, Nagy F, et al. ELONGATED HYPOCOTYL 5 mediates blue light signalling to the Arabidopsis circadian clock. Plant J. 2018;96(6):1242–54. https://doi.org/10.1111/tpj.14106.

    CAS  Article  PubMed  Google Scholar 

  49. 49.

    Osnato M, Matias-Hernandez L, Aguilar-Jaramillo AE, Kater MM, Pelaz S. Genes of the RAV family control heading date and carpel development in Rice. Plant Physiol. 2020;183(4):1663–80. https://doi.org/10.1104/pp.20.00562.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  50. 50.

    Neff MM, Chory J. Genetic interactions between phytochrome a, phytochrome B, and cryptochrome 1 during Arabidopsis development. Plant Physiol. 1998;118(1):27–35. https://doi.org/10.1104/pp.118.1.27.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  51. 51.

    Fang J, Zhang F, Wang H, Wang W, Zhao F, Li Z, et al. Ef-cd locus shortens rice maturity duration without yield penalty. Proc Natl Acad Sci U S A. 2019;116(37):18717–22. https://doi.org/10.1073/pnas.1815030116.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  52. 52.

    Bazile D, Jacobsen SE, Verniau A. The global expansion of quinoa: trends and limits. Front Plant Sci. 2016;7:622.

    PubMed  PubMed Central  Google Scholar 

  53. 53.

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

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  54. 54.

    Ding Z, Tie W, Fu L, Yan Y, Liu G, Yan W, et al. Strand-specific RNA-seq based identification and functional prediction of drought-responsive lncRNAs in cassava. BMC Genomics. 2019;20(1):214. https://doi.org/10.1186/s12864-019-5585-5.

    Article  PubMed  PubMed Central  Google Scholar 

  55. 55.

    Ma X, Zhang X, Traore SM, Xin Z, Ning L, Li K, Zhao K, Li Z, He G, Yin D: Genome-wide identification and analysis of long noncoding RNAs (lncRNAs) during seed development in peanut (Arachis hypogaea L.). BMC Plant Biol 2020, 20(1):192.

    CAS  Article  Google Scholar 

  56. 56.

    Wickland DP, Hanzawa Y. The FLOWERING LOCUS T/TERMINAL FLOWER 1 gene family: functional evolution and molecular mechanisms. Mol Plant. 2015;8(7):983–97. https://doi.org/10.1016/j.molp.2015.01.007.

    CAS  Article  PubMed  Google Scholar 

  57. 57.

    Komiya R, Ikegami A, Tamaki S, Yokoi S, Shimamoto K. Hd3a and RFT1 are essential for flowering in rice. Development. 2008;135(4):767–74. https://doi.org/10.1242/dev.008631.

    CAS  Article  PubMed  Google Scholar 

  58. 58.

    Liu W, Jiang B, Ma L, Zhang S, Zhai H, Xu X, et al. Functional diversification of flowering locus T homologs in soybean: GmFT1a and GmFT2a/5a have opposite roles in controlling flowering and maturation. New Phytol. 2018;217:1335–45. https://doi.org/10.1111/nph.14884.

  59. 59.

    Chab D, Kolar J, Olson MS, Storchova H. Two flowering locus T (FT) homologs in Chenopodium rubrum differ in expression patterns. Planta. 2008;228(6):929–40. https://doi.org/10.1007/s00425-008-0792-3.

    CAS  Article  PubMed  Google Scholar 

  60. 60.

    Golicz AA, Steinfort U, Arya H, Singh MB, Bhalla PL. Analysis of the quinoa genome reveals conservation and divergence of the flowering pathways. Funct Integr Genomics. 2020;20(2):245–58. https://doi.org/10.1007/s10142-019-00711-1.

    CAS  Article  PubMed  Google Scholar 

  61. 61.

    Tsuji H. Taoka K-i, Shimamoto K: Florigen in rice: complex gene network for florigen transcription, florigen activation complex, and multiple functions. Curr Opin Plant Biol. 2013;16(2):228–35. https://doi.org/10.1016/j.pbi.2013.01.005.

    CAS  Article  PubMed  Google Scholar 

  62. 62.

    Higuchi Y, Sumitomo K, Oda A, Shimizu H, Hisamatsu T. Day light quality affects the night-break response in the short-day plant chrysanthemum, suggesting differential phytochrome-mediated regulation of flowering. J Plant Physiol. 2012;169(18):1789–96. https://doi.org/10.1016/j.jplph.2012.07.003.

    CAS  Article  PubMed  Google Scholar 

  63. 63.

    Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35(Web Server issue):W345–9.

    Article  Google Scholar 

  64. 64.

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

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  65. 65.

    Wang L, Park HJ, Dasari S, Wang S, Kocher JP, Li W. CPAT: coding-potential assessment tool using an alignment-free logistic regression model. Nucleic Acids Res. 2013;41(6):e74. https://doi.org/10.1093/nar/gkt006.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  66. 66.

    Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, Mitchell AL, et al. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 2016;44(D1):D279–85. https://doi.org/10.1093/nar/gkv1344.

    CAS  Article  PubMed  Google Scholar 

  67. 67.

    Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30. https://doi.org/10.1093/nar/28.1.27.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  68. 68.

    Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The gene ontology consortium. Nat Genet. 2000;25(1):25–9. https://doi.org/10.1038/75556.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  69. 69.

    Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504. https://doi.org/10.1101/gr.1239303.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

Download references

Funding

This work is supported by the National Natural Science Foundation of China (Grant 31701493); the Sichuan Science and Technology Program (Grant 2018JY0344, Grant 2020YJ0199); the China Agriculture Research System (Project CARS-08-02A). The funding bodies have no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information

Affiliations

Authors

Contributions

QW conceived and designed this study. QW wrote the manuscript. QW, YL and XW performed most of the experiments. XY, CL, XB and QL participated in analyzing the RNA-seq data. YL, DX and YW prepared the samples for sequencing. LZ and GZ participated in quinoa flowering phenotype analysis. All authors have read and approved the final manuscript.

Corresponding author

Correspondence to Qi Wu.

Ethics declarations

Ethics approval and consent to participate

The plant seeds used in this study were harvested from the Experimental Station of Key Laboratory of Coarse Cereal Processing, Ministry of Agriculture and Rural Affairs. The plants were grown in green house and controlled growth chamber, and all the experiments were carried out in Key Laboratory of Coarse Cereal Processing, Ministry of Agriculture and Rural Affairs, Chengdu University. All the Experimental research and field studies, including the collection of plant material, comply with the guidelines of Key Laboratory of Coarse Cereal Processing, Ministry of Agriculture and Rural Affairs, Chengdu University.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Fig. S1.

Spearman correlation coefficients between different samples

Additional file 2: Table S1.

Expression levels of all DE lncRNAs

Additional file 3: Table S2.

Expression levels of all DEGs

Additional file 4: Table S3.

Targets of the up-regulated lncRNAs in SD1_vs_SD2

Additional file 5: Table S4.

Targets of the down-regulated lncRNAs in SD1_vs_SD2

Additional file 6: Table S5.

Targets of the up-regulated lncRNAs in SD2_vs_NB

Additional file 7: Table S6.

Targets of the down-regulated lncRNAs in SD2_vs_NB

Additional file 8: Table S7.

Network of the 17 positive flowering lncRNAs

Additional file 9: Table S8.

Network of the 7 negative flowering lncRNAs

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Wu, Q., Luo, Y., Wu, X. et al. Identification of the specific long-noncoding RNAs involved in night-break mediated flowering retardation in Chenopodium quinoa. BMC Genomics 22, 284 (2021). https://doi.org/10.1186/s12864-021-07605-2

Download citation

Keywords

  • LncRNA
  • Night break
  • Flowering
  • Chenopodium quinoa
  • Co-expression network