- Open Access
Genome-wide profiling of long noncoding RNAs involved in wheat spike development
BMC Genomics volume 22, Article number: 493 (2021)
Long noncoding RNAs (lncRNAs) have been shown to play important roles in the regulation of plant growth and development. Recent transcriptomic analyses have revealed the gene expression profiling in wheat spike development, however, the possible regulatory roles of lncRNAs in wheat spike morphogenesis remain largely unclear.
Here, we analyzed the genome-wide profiling of lncRNAs during wheat spike development at six stages, and identified a total of 8,889 expressed lncRNAs, among which 2,753 were differentially expressed lncRNAs (DE lncRNAs) at various developmental stages. Three hundred fifteen differentially expressed cis- and trans-regulatory lncRNA-mRNA pairs comprised of 205 lncRNAs and 279 genes were predicted, which were found to be mainly involved in the stress responses, transcriptional and enzymatic regulations. Moreover, the 145 DE lncRNAs were predicted as putative precursors or target mimics of miRNAs. Finally, we identified the important lncRNAs that participate in spike development by potentially targeting stress response genes, TF genes or miRNAs.
This study outlines an overall view of lncRNAs and their possible regulatory networks during wheat spike development, which also provides an alternative resource for genetic manipulation of wheat spike architecture and thus yield.
The architecture of flowering plants is largely influenced by the varied inflorescences and their spatial positions and orientations . The inflorescence development begins with the transition of a vegetative shoot apical meristem (SAM) to a reproductive inflorescence meristem (IM) after perceiving intrinsic and extrinsic signals, such as photoperiod and temperature [2, 3]. In the model eudicot plant Arabidopsis, the IM directly generates floral meristem (FM), which subsequently produces floral organ primordia and eventually develops into a simple raceme-type inflorescence . By contrast, in monocot crop such as wheat, the IM initiates spikelet meristem (SM) and FM, and finally forms a raceme-like unbranched inflorescence named spike. Although accumulating evidence has demonstrated that some key regulators governing inflorescence development, including the well-known ABCDE model genes for floral organ development, are functionally conserved among different plant species [5, 6], the largely diversified morphology of inflorescences in plant kingdom suggests that some divergent regulators or mechanisms may exist and contribute to specific inflorescence architecture among different plant species and crops.
As the reproductive organ, wheat spike is tightly related to agronomic traits and thus yield. For example, the grain number per spike is one of the crucial determinants of the yield . The wheat spike development involves the subsequent formation and development of SM, glume primordium, FM, stamen and pistil primordia, and floral organs . Recently, several important regulators governing wheat spike development have been characterized. For examples, the TaTFL1-2D, an ortholog of Arabidopsis TERMINAL FLOWER 1 (TFL1), regulates wheat spike complexity by altering the number of spikelet, floret, and grain per spike , while the VERNALIZATION 1 (VRN1, also named as FRUITFULL 1 (FUL1), FUL2 and FUL3 function redundantly in determining SM identity, as the lateral meristems in the vrn1 ful2 ful3 triple mutant are reverted to vegetative meristems to form leaves instead of spikelets . The FRIZZY PANICLE (FZP) modifies the spike structure by regulating SM identity, since its mutation results in the ectopic formation of SM and thus supernumerary spikelets . Moreover, recent studies have begun to reveal the gene expression profiling during wheat spike development. The transcriptome of wheat inflorescence highlighted the gene regulatory networks between miRNAs and their target genes , and the gene expression profiles at six developmental stages of wheat spike revealed the expression dynamics of some important genes regulating spike development, including VRN1, Photoperiod-1 (Ppd-1) and Earliness per se 3 (Eps-3) . An associative transcriptome analysis of 90 wheat varieties also identified a few of genes related to spike complexity .
LncRNAs are defined as the RNA transcripts more than 200 base pairs that lack protein-encoding potential, which are less abundant but exhibit more spatiotemporally specific expression patterns than the protein-coding genes . Previous studies have shown that lncRNAs can function in a cis manner by affecting neighboring loci, and in trans by performing distal regulatory functions, as scaffolds for protein complexes, decoys, guides, and enhancers to regulate gene expression [15, 16]. Moreover, some of lncRNAs are reported to function as the precursors of miRNAs to influence miRNA biosynthesis or as mimic targets and thus interfering miRNA functions . In Arabidopsis and rice, lncRNAs have been reported to play essential roles in the development of reproductive organs, including floral transition, meiosis progression, anther and pollen development [18,19,20,21,22]. For examples, the Arabidopsis COOLAIR and COLDAIR influence flowering time by regulating FLOWERING LOCUS C (FLC) expression at epigenetic and posttranscriptional levels , while the asHSFB2a affects female gametophyte development . The rice long-day-specific male-fertility-associated RNA (LDMAR) regulates anther development and male sterility in a photoperiod-sensitive manner . Recent works also suggest that lncRNAs participate in grain development and stress response in wheat [24, 25]. However, the regulatory roles of lncRNAs during wheat spike development remain largely elusive.
Here, we analyzed the dynamic transcriptome of developing wheat spike, and identified a total of 8,889 expressed lncRNAs with 2,753 DE lncRNAs at different developmental stages. We further predicted 315 DE lncRNA-mRNA regulatory pairs and analyzed the molecular events regulated by these DE lncRNA. Subsequently, 24 important lncRNAs and their potential regulatory genes were identified. We also predicted some lncRNAs may function as putative precursors or target mimics of miRNAs. Our data thus provide an overall view on dynamic expression and potential regulation of lncRNAs during wheat spike development, which are also alternative resources for further manipulation of wheat spike architecture and thus yield.
The lncRNAs transcriptomic profiling of developing wheat spike
To monitor the dynamic profiling of lncRNAs during wheat spike development, the developing spikes of winter wheat Zhengmai366 were collected at six stages (S1-S6) referred as IM stage (S1), SM stage (S2), glume primordium stage (S3), FM stage (S4), stamen and pistil primordium stage (S5) and floral organ stage (S6), respectively, according to anatomic and morphological features described previously  (Supplementary Fig. S1). These developmental stages morphologically corresponded to the W2.0 (S1), W2.5 (S2), W3.0 (S3), W3.5 (S4), W5 (S5) and W6 (S6) on the respective Waddington scale as described previously .
A total of 18 RNA samples from six developmental stages with three biological replicates were subjected to RNA isolation and sequencing, respectively. Approximately, 15 Gb clean base pairs for each sample were aligned to the wheat genomic sequence (IWGSC Refseq v1.0, Ensemble Plants; Supplementary Table S1). Unique reads were used to quantify gene models in fragments per kilobase of exon model per million mapped reads (FPKM). With average FPKM ≥ 1 of three replicates in at least one stage, a total of 48,345 coding genes were expressed, among which 11,792 were differentially expressed during spike development (Supplementary Table S2; Supplementary Table S3), consistent with ~ 50,000 genes detected in previous study .
Following LGC (https://bigd.big.ac.cn/lgc/), Coding-Non-Coding-Index (CNCI) predictions and Pfam blast of potential coding capacity, a total of 8,889 transcripts with average FPKM ≥ 1 in at least one stage were considered as the high confidence lncRNAs expressed during spike development, among which 143 had been annotated in the total of 362 lncRNAs annotated in wheat genome (Supplementary Table S4). For individual developmental stage, there were ~ 8,000 lncRNAs, including ~ 3,000 lncRNAs with low abundance (transcripts per kilobase of exon model per million mapped reads, TPM value between 1 ~ 10), ~ 4,500 lncRNAs with medium abundance (TPM value between 10 ~ 100), and ~ 500 lncRNAs with high abundance (TPM value > 100) (Fig. 1a). These expressed lncRNAs were almost evenly distributed across all chromosomes of the wheat genome (Fig. 1b). However, the subgenome distribution showed that lncRNAs appeared to be slightly preferential on B subgenome (A, 2,677, 31.52 %; B, 3,308, 38.95 %; D, 2,509, 29.54 %; Fig. 1c). Moreover, the lncRNAs expression levels between any two developmental stages were highly correlated with the correlation coefficient above 0.98 (Fig. 1d), consistent with the continuity of spike development despite of distinct morphological characteristics existed at each stage. Importantly, among the expressed lncRNAs, 7,522 (84.7 %) were detected in all the six developmental stages (Supplementary Fig. S2). The numbers of lncRNAs expressed at only one stage (FPKM < 1 at the other stages) were 63, 70, 60, 43, 41 and 61 from S1 to S6, respectively (Supplementary Fig. S2; Supplementary Table S5). These observations demonstrate that lncRNAs are highly involved in spike development and some of them may function in a stage-specific manner to contribute to specific developmental characteristic.
Expression clustering of differentially expressed lncRNAs
To explore the possible regulation of lncRNAs during spike development, we next focused on the lncRNAs that were differentially expressed among the six developmental stages. Through pair-wise differential analysis by DEseq2 , a total of 2,753 DE lncRNAs with at least two-fold change in expression and a p-value less than 0.05 were identified (Supplementary Table S6). The numbers of DE lncRNAs between the earlier stages (S1-S3) were relatively less than those between the earlier stage (S1-S3) and later stage (S5-S6) (Supplementary Fig. S3), consistent with the progressive complexity of spike architecture during development. Using k-means clustering, we divided all DE lncRNAs into nine expression clusters (Fig. 2; Supplementary Table S6). The 324 lncRNAs in cluster 1 were those whose expressions were consistently elevated from S1 to S6, while 655 lncRNAs in clusters 2 and 3 included those whose expressions were consistently decreased from S1 to S6 with varied expression trends (Fig. 2). These clusters largely represented the lncRNAs that might function as positive or negative regulators during whole spike developmental progression. In contrast, the clusters 4-9 represented the lncRNAs that were abundantly expressed at a specific developmental stage among S1-S6 (Fig. 2). These lncRNAs likely participate in regulation of the distinct molecular events occurred at specific developmental stage. Further qRT-PCR analysis of six selected lncRNAs predominantly expressed at each developmental stage, including MSTRG.60752, MSTRG.38633, MSTRG.2756, MSTRG.34036, MSTRG.11486, and MSTRG.39458, validated the high reliability of their expression trends revealed by transcriptome data (Supplementary Fig. S4). In addition, we also amplified four lncRNAs, MSTRG.59353, MSTRG.60752, MSTRG.38633 and MSTRG.11486, and validated their identity and expression through sanger sequencing.
The coding genes potentially targeted by DE lncRNAs
Since lncRNAs can regulate gene expression either in cis or trans manner, we next analyzed the differentially expressed coding genes (DEGs) potentially targeted by DE lncRNAs based on the position and sequence relationships as well as their expression correlation coefficients. The lncRNAs and their potential target coding genes were thus referred as DE lncRNA-mRNA pairs. A total of 315 DE lncRNA-mRNA pairs composed of 205 lncRNAs and 279 coding genes were identified, of which 204 pairs were predicted to be cis-regulatory and 111 were trans-regulatory (Fig. 3a; Supplementary Table S7). Interestingly, all these lncRNA-mRNA pairs had a positive correlation in their expression patterns (Supplementary Table S7). Among these regulatory pairs, more than 80 % of lncRNAs (166) targeted the individual coding gene, while 29 lncRNAs targeted two or three genes, and ten lncRNAs had more than five target genes (Fig. 3b). Likewise, more than 92 % of coding genes (257) corresponded to an individual lncRNA, while 20 genes were targeted by two or three lncRNAs, and only two genes were targeted by four and five lncRNAs (Fig. 3c). These observations implicate that the regulatory crosstalk occurs between lncRNAs and coding genes. Further qRT-PCR analyses of three representative lncRNAs, MSTRG.21315, MSTRG.42234, MSTRG.28564 and their predicted targets verified their expression correlations during spike development (Fig. 3d-f).
GO enrichment analysis of coding genes regulated by DE lncRNAs
To explore the molecular events regulated by DE lncRNAs, we performed Gene Ontology (GO) analysis of their target coding genes. Interestingly, among all the DEGs potentially targeted by DE lncRNAs, the genes associated with “response to stress” and “response to chemical” were highly enriched in biological process, and the genes associated with “transcription regulator activity” and “enzyme activity” were enriched in molecular function (Fig. 4a). Consistent with this, GO enrichment analysis of coding genes targeted by lncRNAs with gradually increased expression during spike development in cluster 1 further revealed that genes related to “response to wounding” and “nucleus” were enriched (Fig. 4b), and that genes associated with “macromolecule methylation” and “homeostatic process” were enriched among the DEGs targeted by lncRNAs with gradually decreased expression in cluster 2 and 3 (Fig. 4c). Therefore, the regulatory roles of lncRNAs during wheat spike development seemed to be mainly associated with stress response, transcriptional and enzymatic regulation.
As the coding genes targeted by stage-specific lncRNAs might contribute to the spike patterning at an individual developmental stage, we also examined the genes targeted by lncRNAs in clusters 4–9 (Fig. 5; Supplementary Table S8). As expected, these coding genes indeed exhibited predominant expressions at the corresponding developmental stage as did their paired lncRNAs (Fig. 5). At S1, there were four genes encoding metabolic related proteins, three encoding Ycf68 proteins with unknown functions, several genes encoding cytochrome, glutamate receptor, pentatricopeptide repeat-containing protein and zinc finger protein (Fig. 5a). At S2, three lncRNAs were involved in regulation of three genes encoding a trigger factor, a ribosomal protein, and an auxin response factor (Fig. 5b). At S3, the genes encoding chymotrypsin inhibitors, and sugar or lipid metabolism related proteins were identified as potential targets of lncRNAs (Fig. 5c). Interestingly, a large portion of genes encoding wound-responsive proteins, dehydrin proteins and early nodulin related proteins were found to be targeted by lncRNAs at S4 (Fig. 5d). Similarly, among 33 genes targeted by lncRNAs at S5, there were 22 genes that encode wound-responsive proteins (Fig. 5e; Supplementary Table S8). This observation implicates that lncRNAs-regulated stress response especially wounding response may contribute to floret, stamen and pistil primordium formation at S4 and S5. Among genes potentially targeted by lncRNAs at S6, there were several genes encoding MADS-box TFs, zinc finger proteins and Argonaute1 proteins, and a large portion of genes were found to encode nuclear proteins, such as chromatin related proteins (Fig. 5f), implying that the transcription regulation and epigenetic modification regulated by lncRNAs are critical for floral organ development.
Identification of important lncRNAs involved in wheat spike development
To identify the important lncRNAs involved in spike development, we closely examined the DE lncRNAs related with stress response, TFs regulation, hormone signaling and chromatin modeling (Table 1; Fig. 5). Among DE lncRNAs related to stress response, there were six lncRNAs, including MSTRG.18189, MSTRG.18191, MSTRG.14917, MSTRG.11521, MSTRG.14912 and MSTRG.11500, that potentially targeted a total of 30 wound-responsive genes (Table 1), while MSTRG.59826 and MSTRG.53898 were found to target nine nodulin-related genes, and MSTRG.47943 targeting five dehydrin genes (Table 1; Fig. 5). By contrast, there were 12 TF genes identified in the DE lncRNA-mRNA pairs, including two MADS-box, eight bZIP and two MYB genes (Table 1). Noticeably, MSTRG.28564 and MSTRG.61765 potentially targeted two MADS-box genes, TraesCS4A02G078700 and TraesCS7D02G380300, while MSTRG.51463, MSTRG.62263, MSTRG.39593 and MSTRG.43002 potentially targeted eight bZIP genes, and MSTRG.57306, MSTRG.30797 targeted two MYB genes (Table 1). Moreover, auxin has been reported to play critical roles during inflorescence development , we also observed that four auxin response genes, including TraesCS2D02G491200, a homolog of AUXIN RESPONSE FACTOR 5 (ARF5) in Arabidopsis, were targeted by four lncRNAs (Table 1). In addition, Argonaute1 (AGO1) can bind to chromatin to regulate transcription in association with SWI/SNF complexes . We also found that MSTRG.59353 could target AGO1d-7A and AGO1d-7B, and that MSTRG.23144 and MSTRG.19854 possibly targeted two SWI/SNF-related genes (Table 1). These DE lncRNAs and their target genes likely represent the key regulatory networks of lncRNAs during wheat spike patterning and development.
LncRNAs as putative miRNA precursors or target mimics
LncRNAs can also act as precursors or target mimics of miRNA to regulate miRNAs’ function, and some of miRNAs have been reported to be involved in wheat spike development, such as miR172 and miR156 [29,30,31]. To identify the lncRNAs that possibly act as precursors or target mimics of miRNA, we further performed sequence alignment to miRBase and psRNAtarget [32, 33]. A total of 58 DE lncRNAs were predicted as putative precursors of 21 miRNAs, including tae-miR156, tae-miR160, tae-miR167 and tae-miR444 (Supplementary Table S9). We also observed that one lncRNA might serve as the precursors for one or a few of miRNAs, while several lncRNAs might be potential precursors for the same miRNA (Supplementary Table S9). Notably, MSTRG.41996 and MSTRG.45365 that were highly expressed at S1 and S2 were predicted as the precursors of tae-miR167a and tae-miR167b, while MSTRG.41907 highly expressed during S1-S3 was predicted as the precursor of tae-miR160, and MSTRG.50297 predominantly expressed at S1 and S6 was predicted as the precursor of tae-miR444a (Supplementary Table S6, S9). To further verify this, we used RNAfold web server (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi) to predict the secondary structures of the two representative lncRNAs and their miRNA precursors. As shown in Fig. 6a, the secondary structure of MSTRG.41907.1 exhibited multiple stem-loop structures, one of which could be cleaved to release the precursor sequence of tae-miR160 and eventually form mature tae-miR160. Similarly, the secondary structure of MSTRG.50297.1 might be cleaved to produce precursor of tae-miR444a and form mature tae-miR444a (Fig. 6b). Furthermore, transcripts from 87 DE lncRNAs were predicted as target mimics of 37 miRNAs, including tae-miR1139, tae-miR9672, tae-miR9674, tae-miR9783 and other miRNAs (Supplementary Table S10). These observations support that these lncRNAs participate in wheat spike development as precursors or target mimics of miRNA.
LncRNAs have been reported to play important roles in plant growth and development in various plant species [18, 20,21,22]. Recent works have also identified a series of lncRNAs during wheat grain development, tillering development, and defense-response [25, 34, 35]. However, lncRNAs and their possible regulations during wheat spike development remain elusive. Here, we monitored the genome-wide profiling of lncRNAs and identified a total of 8,889 lncRNAs expressed in wheat spike, of which 2,753 lncRNAs are differentially expressed among six developmental stages. We predicted 315 DE lncRNA-mRNA regulatory pairs and 145 potential miRNA precursors or target mimics, and outlined the possible molecular events regulated by lncRNAs (Fig. 7). Our work provides a glimpse of lncRNAs and their molecular regulations during wheat SM, FM and floral organ development. As wheat spike development is critical for determination of agronomic traits such as spikelet, floret and grain numbers, further works on these lncRNAs and their potential targets will be necessary not only for understanding the lncRNAs and their molecular regulations during wheat spike development, but also for genetic modification of spike architecture and thus wheat yield.
Since lncRNAs can execute their functions by cis- or trans-regulation of target coding genes [15, 16], the identification of their possible target genes is necessary for understanding their regulatory roles. Interestingly, our GO enrichment of genes targeted by DE lncRNAs revealed that these genes were mainly associated with stress response, transcription regulation and enzyme activity regulation (Fig. 4), suggesting that lncRNAs may mainly target on these biological processes to regulate spike development. Indeed, stress responses in plants, such as drought, cold, hormone and so on, not only influence wheat spikelet and grain numbers and thus yield under ever-changed environments [36,37,38], but also play important roles in the cell fate determination during plant meristem differentiation and regeneration [39, 40]. For example, wounding is one of key signals to trigger cell reprograming in plant regeneration, and the following auxin, jasmonic acid, ethylene, and electrical pulses are critical for wound-triggered cell fate switches [41, 42]. As wheat spike formation and development involve multiple meristem differentiation and organ formation, it is likely that lncRNAs may participate in regulation of such meristem differentiation by targeting the stress and hormone related genes. Consistent with this, four genes encoding auxin response factors were regulated by four lncRNAs with various expression patterns (Table 1). Meanwhile, a total of 134 and 27 DEGs related with auxin and cytokinin were also detected during spike development, respectively (Supplementary Table S3). On the other hand, TFs, such as AP2, SPL, bZIP and MADS-box, have been reported to play essential roles during inflorescence development in wheat and other plants [5, 6, 10, 30]. Here, we found that 12 TF genes were potentially targeted by lncRNAs (Table 1), suggesting that these lncRNAs may also facilitate their roles on transcriptional regulation of key TFs required for wheat spike development. In addition, metabolic process influences spike and grain numbers in crops mainly through accumulation of assimilates, especially for carbohydrates, which is controlled by a series of enzymes [43,44,45,46]. The enrichment of enzyme-related genes demonstrates that cellular metabolic processes essential for spike organs (spikelet, glume and floret) development are also regulated by lncRNAs.
As wheat spike development undergoes subsequent developmental processes of multiple meristems and floral organs (S1-S6), identification of the stage-specific lncRNAs and their target genes will be helpful to understand the specific events for spike patterning regulated by lncRNAs. We observed that some of metabolic related genes targeted by lncRNAs were mainly enriched during S1-S4 (Fig. 5; Supplementary Table S8), suggesting that the regulation of these metabolic processes by lncRNAs may contribute to the formation of various meristems (spikelet, glume, lemma, and floret) at early developmental stages. By contrast, the enrichment of a series of genes associated with wound responses by lncRNAs at S4 and S5 implicates that the floret, stamen and pistil primordium formation may require the precise regulations of lncRNAs during the cell fate changes and environmental responses at these stages (Fig. 5; Supplementary Table S8). Moreover, the enrichment of TF and chromatin-related genes targeted by lncRNAs at S5 and S6 demonstrates that lncRNAs may regulate floral organ development through transcriptional and epigenetic modification (Fig. 5; Supplementary Table S8). Finally, we identified 24 important lncRNAs that were involved in various regulations during wheat spike development and found some of them may mainly function at specific stage (Table 1). For examples, MSTRG.18189, MSTRG.11521, MSTRG.14912 and MSTRG.14917 could target a series of wound responsive genes at S4 and S5, and MSTRG.59353 possibly target AGO1d-7 A and AGO1d-7B at S6 (Fig. 5; Table 1). As wheat AGO1d has been shown to play an important role in anther and pollen development , it is likely that MSTRG.59353 may be involved in floret fertility development. Therefore, these important lncRNAs and their potential targets will be the candidates for further elucidating the molecular regulations of lncRNAs and for genetically modification of wheat spike architecture.
MiRNAs play important roles in regulating plant morphogenesis, including the transition of reproductive phase from vegetative phase and the establishment of organ patterning [47,48,49,50]. In wheat, several key miRNAs have been reported to participate in spike development. For examples, miRNA156 targets the SPL gene to regulate spikelet development, while miRNA172 modifies wheat spike architecture through Q gene [31, 51, 52]. Transcriptomic analysis also found that miR159, miR167, and miR319 were abundantly expressed in wheat inflorescence meristems . Here, we also predicted some DE lncRNAs as putative precursors or target mimics of these important miRNAs, strengthening that these lncRNAs may facilitate their regulations in wheat spike development through modifying functional miRNAs and their target genes. For examples, MSTRG.2923 may act as the precursor of tae-miR156 and tae-miR160, and MSTRG.45365 as the precursor of tae-miR167 (Supplementary Table S9). We also identified some lncRNAs may function as precursors or targets of miRNAs with unknown functions, such as miR1135 and miR1136. Further works on these lncRNAs and their corresponding miRNAs will be helpful to clarify their molecular interactions during wheat spike development.
In this study, a total of 8,889 lncRNAs are identified to be expressed during wheat spike development, among which 2,753 DE lncRNAs are categorized into nine expression clusters. A total of 315 DE lncRNA-mRNA cis- and trans-regulatory pairs are predicted, among which 24 important lncRNAs and their regulatory genes are identified. They are mainly related to stress responses, transcriptional regulation, and enzyme activity regulation. Moreover, 58 DE lncRNAs are identified as potential miRNA precursors and 87 DE lncRNAs as miRNA target mimics. These findings provide an overall view on lncRNAs and their possible regulatory networks in wheat spike development.
Plant materials and sample collection
The winter wheat cultivar Zhengmai366, a widely bred cultivar for high-yield and good-quality wheat, is obtained from the breeder of Wheat Research Institute, Henan Academy of Agricultural Sciences, with permissions for usage and collection of plants and seeds. All the plant growth management and methods used in this study are complied with relevant national guidelines and legislations of China (http://www.moa.gov.cn). The Zhengmai366 was grown in Beijing, China and the developmental stages of spikes were determined based on the anatomic and morphological features under stereomicroscope OLYMPUS SZX16 (Japan) according to previous study . The developing spikes at six developmental stages referred as S1 to S6 were collected, which represented the developmental stages of IM, SM, glume primordium, FM, stamen and pistil primordia, and floral organs (anther and awn), respectively. The dissected spikes were immediately frozen in liquid nitrogen and stored at -80 °C. About 20-50 spikes at each developmental stage were pooled for each of the three biological replicates and subjected for RNA isolation and sequencing.
Scanning electron microscopy
To determine the developmental stage, the developing spikes were fixed in FAA (50 % v/v ethanol, 5 % v/v acetic acid, and 3.7 % v/v formaldehyde) for 24 h at 4 °C, dehydrated through a standard ethanol series, and dried with CO2 critical point. The dried samples were mounted on stubs, coated with gold, and photographed with a Hitachi S-4800 device (Tokyo, Japan).
RNA isolation, library construction and sequencing
Total RNA was isolated with TRIzol reagent (Invitrogen) and treated with DNaseI (Invitrogen) to eliminate contaminating DNA. Construction of transcriptome libraries and deep sequencing were performed by The Beijing Genomics Institute. Ribosomal RNA was removed using a Ribo-Zero™ rRNA Removal Kit (Illumina). Subsequently, random primers from TruSeq® Stranded Kit (Illumina) were used to synthesize cDNA with the templates of fragmented RNAs. The cDNA library was obtained by PCR amplification. The resulting libraries were sequenced on a HiSeq 4000 instrument (Illumina) that generated paired end reads of 100 bp.
Reads mapping, transcript assembly and lncRNAs identification
After removing the adaptor-polluted and low-quality reads, approximately 15 Gb clean base pairs for each sample were mapped to the Triticum aestivum L. reference genome assembly (IWGSC Refseq v1.0, Ensemble Plants, http://ftp.ensemblgenomes.org/pub/plants/release-46/) using HISAT with default parameters . Stringtie was used to assemble the transcripts and Gffcompare was used to compare these transcripts with known gene annotation models to predict novel transcripts . The protein-coding capacity of novel transcripts with length of > 200 bp was identified using LGC and CNCI with default parameters, respectively . Then Pfam database was used to ensure the predicted transcripts have no protein-coding domains . The transcripts as lncRNAs were determined by consistency of the three judgement methods.
Expression analysis of lncRNAs
FPKM were calculated to estimate the expression level of genes in each sample. The lncRNAs and coding genes with an average FPKM of ≥ 1 in at least one stage were considered as expressed. TPM were calculated to compare the gene abundance levels in the same stage. Circos software was used to visually describe the distribution of expressed lncRNAs on 21 chromosomes . DEseq2 was used to perform differential expression analysis among various developmental stages . DE lncRNAs and DEGs were filtered with at least two-fold expression change and a p-value less than 0.05. Clustering of DE lncRNAs was performed by R version 3.6.2 using k-means function. Heat maps and the Venn diagram were drawn using TBtools .
Validation of selected lncRNAs by qRT-PCR and Sanger sequencing
Total RNA was isolated with TRIzol reagent (Invitrogen). Reverse transcription was performed using the PrimeScript RT reagent Kit with gDNA Eraser following the manufacturer’s recommendations (Takara). Real-time quantitative PCR was performed using the SYBR Premix Ex TaqTM Kit (Takara) on a LightCycler 96 machine (Roche). The wheat ACTIN gene was used as a reference control and relative expression levels were calculated using the 2-▵▵CT method . The qRT-PCR was performed with three biological replicates and all the primers sequence information was listed in Supplementary Table S11. Four lncRNAs were amplified using the Phusion high-fidelity enzyme (NEB) and sequenced by Sanger sequencing in The Beijing Genomics Institute.
Prediction of lncRNA target genes and GO analysis
The lncRNA and its predicted cis- or trans-regulatory coding genes were considered as lncRNA-mRNA pairs. The cis-regulatory genes of lncRNA were determined if there was an overlap between the lncRNA and the coding genes, or the lncRNA located within 100 kb up- or down-stream of the coding genes . For trans-regulation, LncTar was used to calculate the free energy of lncRNA-mRNA pairs and the normalized free energy at -0.1 was set as cutoff to determine target genes . In addition, the lncRNA-mRNA pairs were filtered with the Pearson correlation efficient of their expression levels (|r|> 0.9). The DEGs possibly targeted by DE lncRNAs were considered as candidate coding genes of interest for further analysis. GO enrichment analysis was performed with TBtools software. The GO terms with a p-value below 0.05 were considered significantly enriched.
Prediction of miRNA precursors or target mimics
To identify potential miRNA precursors, lncRNAs were blasted to miRBase (http://www.mirbase.org/search.shtml) and the identity between lncRNAs and miRNA precursors that were 100 % and e-value < 0.001 were selected . The miRNA target mimics prediction was performed by aligning the wheat mature miRNA sequences against lncRNA sequences using psRNAtarget with default parameters except for a strict Expectation value 2 .
Availability of data and materials
The raw sequence reads are available in NCBI sequence read archive database and can be accessed under the ID PRJNA718695 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA718695) (NCBI Sequence Read Archive SRR14116504, SRR14116505, SRR14116506, SRR14116507, SRR14116508, SRR14116509, SRR14116510, SRR14116511, SRR14116512, SRR14116513, SRR14116514, SRR14116515, SRR14116516, SRR14116517, SRR14116518, SRR14116519, SRR14116520 and SRR14116521). All other relevant data generated or analyzed in this study are included in this article and its Additional files.
Differentially expressed coding genes
- DE lncRNAs:
Differentially expressed lncRNA
Fragments per kilobase of exon model per million mapped reads
Long noncoding RNAs
Shoot apical meristem
Transcripts per kilobase of exon model per million mapped reads
Vollbrecht E, Springer PS, Goh L, Buckler ESt, Martienssen R. Architecture of floral branch systems in maize and related grasses. Nature. 2005;436(7054):1119–1126.
Boss PK, Bastow RM, Mylne JS, Dean C. Multiple pathways in the decision to flower: enabling, promoting, and resetting. Plant Cell. 2004;16(Suppl):S18-31.
Wigge PA, Kim MC, Jaeger KE, Busch W, Schmid M, Lohmann JU, Weigel D. Integration of spatial and temporal information during floral induction in Arabidopsis. Science. 2005;309(5737):1056–1059.
Smyth DR, Bowman JL, Meyerowitz EM. Early flower development in Arabidopsis. Plant Cell. 1990;2(8):755–767.
Ali Z, Raza Q, Atif RM, Aslam U, Ajmal M, Chung G. Genetic and molecular control of floral organ identity in cereals. Int J Mol Sci. 2019;20(11):2743.
Callens C, Tucker MR, Zhang D, Wilson ZA. Dissecting the role of MADS-box genes in monocot floral development and diversity. J Exp Bot. 2018;69(10):2435–2459.
Sreenivasulu N, Schnurbusch T. A genetic playground for enhancing grain number in cereals. Trends Plant Sci. 2012;17(2):91–101.
Bonnett OT. The development of the wheat spike. J Agric Res. 1936;53:445–451.
Wang Y, Yu H, Tian C, Sajjad M, Gao C, Tong Y, Wang X, Jiao Y. Transcriptome association identifies regulators of wheat spike architecture. Plant Physiol. 2017;175(2):746–757.
Li C, Lin H, Chen A, Lau M, Jernstedt J, Dubcovsky J. Wheat VRN1, FUL2 and FUL3 play critical and redundant roles in spikelet development and spike determinacy. Development. 2019;146(14):dev175398.
Dobrovolskaya O, Pont C, Sibout R, Martinek P, Badaeva E, Murat F, Chosson A, Watanabe N, Prat E, Gautier N et al. FRIZZY PANICLE drives supernumerary spikelets in bread wheat. Plant Physiol. 2015;167(1):189–199.
Feng N, Song G, Guan J, Chen K, Jia M, Huang D, Wu J, Zhang L, Kong X, Geng S et al. Transcriptome profiling of wheat inflorescence development from spikelet initiation to floral patterning identified stage-specific regulatory genes. Plant Physiol. 2017;174(3):1779–1794.
Li Y, Fu X, Zhao M, Zhang W, Li B, An D, Li J, Zhang A, Liu R, Liu X. A genome-wide view of transcriptome dynamics during early spike development in bread wheat. Sci Rep. 2018;8(1):15338.
Li L, Eichten SR, Shimizu R, Petsch K, Yeh CT, Wu W, Chettoor AM, Givan SA, Cole RA, Fowler JE et al. Genome-wide discovery and characterization of maize long non-coding RNAs. Genome Biol. 2014;15(2):R40.
Chekanova JA. Long non-coding RNAs and their functions in plants. Curr Opin Plant Biol. 2015;27:207–216.
Rinn JL, Chang HY. Genome regulation by long noncoding RNAs. Annu Rev Biochem. 2012;81:145–166.
Franco-Zorrilla JM, Valli A, Todesco M, Mateos I, Puga MI, Rubio-Somoza I, Leyva A, Weigel D, Garcia JA, Paz-Ares J. Target mimicry provides a new mechanism for regulation of microRNA activity. Nat Genet. 2007;39(8):1033–1037.
Heo JB, Sung S. Vernalization-mediated epigenetic silencing by a long intronic noncoding RNA. Science. 2011;331(6013):76–79.
Csorba T, Questa JI, Sun Q, Dean C. Antisense COOLAIR mediates the coordinated switching of chromatin states at FLC during vernalization. Proc Natl Acad Sci U S A. 2014;111(45):16160–16165.
Ma J, Yan B, Qu Y, Qin F, Yang Y, Hao X, Yu J, Zhao Q, Zhu D, Ao G. Zm401, a short-open reading-frame mRNA or noncoding RNA, is essential for tapetum and microspore development and can regulate the floret formation in maize. J Cell Biochem. 2008;105(1):136–146.
Wang M, Wu H-J, Fang J, Chu C, Wang X-J. A long noncoding RNA involved in rice reproductive development by negatively regulating osa-miR160. Science Bulletin. 2017;62(7):470–475.
Ding J, Lu Q, Ouyang Y, Mao H, Zhang P, Yao J, Xu C, Li X, Xiao J, Zhang Q. 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–2659.
Wunderlich M, Gross-Hardt R, Schoffl F. Heat shock factor HSFB2a involved in gametophyte development of Arabidopsis thaliana and its expression is controlled by a heat-inducible long non-coding antisense RNA. Plant Mol Biol. 2014;85(6):541–550.
Zhang H, Hu W, Hao J, Lv S, Wang C, Tong W, Wang Y, Wang Y, Liu X, Ji W. Genome-wide identification and functional prediction of novel and fungi-responsive lincRNAs in Triticum aestivum. BMC Genomics. 2016;17:238.
Madhawan A, Sharma A, Bhandawat A, Rahim MS, Kumar P, Mishra A, Parveen A, Sharma H, Verma SK, Roy J. Identification and characterization of long non-coding RNAs regulating resistant starch biosynthesis in bread wheat (Triticum aestivum L.). Genomics. 2020;112(5):3065–3074.
Waddington SR, Cartwright PM, Wall PC. A quantitative scale of spike initial and pistil development in barley and wheat. Ann Bot. 1983;51(1):119–30.
Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26(1):136–138.
Liu C, Xin Y, Xu L, Cai Z, Xue Y, Liu Y, Xie D, Liu Y, Qi Y. Arabidopsis ARGONAUTE 1 binds chromatin to promote gene transcription in response to hormones and stresses. Devl Cell. 2018;44(3):348–361.e347.
Xing S, Salinas M, Hohmann S, Berndtgen R, Huijser P. miR156-targeted and nontargeted SBP-box transcription factors act in concert to secure male fertility in Arabidopsis. Plant Cell. 2010;22(12):3935–3950.
Liu J, Cheng X, Liu P, Sun J. miR156-Targeted SBP-Box transcription factors interact with DWARF53 to regulate TEOSINTE BRANCHED1 and BARREN STALK1 expression in bread wheat. Plant Physiol. 2017;174(3):1931–48.
Debernardi JM, Lin H, Chuck G, Faris JD, Dubcovsky J. microRNA172 plays a crucial role in wheat spike morphogenesis and grain threshability. Development. 2017;144(11):1966–1975.
Griffiths-Jones S, Grocock RJ, van Dongen S, Bateman A, Enright AJ. miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 2006;34(Database issue):D140-144.
Dai X, Zhuang Z, Zhao PX. psRNATarget: a plant small RNA target analysis server (2017 release). Nucleic Acids Res. 2018;46(W1):W49-w54.
Zhang H, Hu W, Hao J, Lv S, Wang C, Tong W, Wang Y, Wang Y, Liu X, Ji W. Genome-wide identification and functional prediction of novel and fungi-responsive lincRNAs in Triticum aestivum. BMC genomics. 2016;17:238–238.
Zhou W, Shi H, Wang Z, Zhao Y, Gou X, Li C, Chen G, Liu S, Deng M, Ma J et al. Identification of lncRNAs involved in wheat tillering development in two pairs of near-isogenic lines. Func Integr Genomic. 2020;20(5):669–679.
Li X, Zhang X, Liu G, Tang Y, Zhou C, Zhang L, Lv J. The spike plays important roles in the drought tolerance as compared to the flag leaf through the phenylpropanoid pathway in wheat. Plant Physiol Bioch. 2020;152:100–111.
Zhang Z, Huang J, Gao Y, Liu Y, Li J, Zhou X, Yao C, Wang Z, Sun Z, Zhang Y. Suppressed ABA signal transduction in the spike promotes sucrose use in the stem and reduces grain number in wheat under water stress. J Exp Bot. 2020;71(22):7241–7256.
Tang Z, Zhang L, Xu C, Yuan S, Zhang F, Zheng Y, Zhao C. Uncovering small RNA-mediated responses to cold stress in a wheat thermosensitive genic male-sterile line by deep sequencing. Plant Physiol. 2012;159(2):721–738.
Xu L. De novo root regeneration from leaf explants: wounding, auxin, and cell fate transition. Curr Opin Plant Biol. 2018;41:39–45.
Ikeuchi M, Favero DS, Sakamoto Y, Iwase A, Coleman D, Rymen B, Sugimoto K. Molecular mechanisms of plant regeneration. Annu Rev Plant Biol. 2019;70:377–406.
Chen L, Sun B, Xu L, Liu W. Wound signaling: the missing link in plant regeneration. Plant Signal Behav. 2016;11(10):e1238548.
Lup SD, Tian X, Xu J, Pérez-Pérez JM. Wound signaling of regenerative cell reprogramming. Plant Sci. 2016;250:178–187.
Fischer RA. Number of kernels in wheat crops and the influence of solar radiation and temperature. J Agr Sci. 1985;105(2):447–461.
Xue GP, McIntyre CL, Jenkins CL, Glassop D, van Herwaarden AF, Shorter R. Molecular dissection of variation in carbohydrate metabolism related to water-soluble carbohydrate accumulation in stems of wheat. Plant Physiol. 2008;146(2):441–454.
Luquet D, Vidal A, Smith M, Dauzat J. ‘More crop per drop’: how to make it acceptable for farmers? Agr Water Manage. 2005;76(2):108–119.
Marcelis LF. Sink strength as a determinant of dry matter partitioning in the whole plant. J Exp Bot. 1996;47 Spec No:1281–1291.
Zhang Z, Zhang X. Argonautes compete for miR165/166 to regulate shoot apical meristem development. Curr Opin Plant Biol. 2012;15(6):652–658.
Spanudakis E, Jackson S. The role of microRNAs in the control of flowering time. J Exp Bot. 2014;65(2):365–380.
Hong Y, Jackson S. Floral induction and flower formation-the role and potential applications of miRNAs. Plant Biotechnol J. 2015;13(3):282–292.
Smoczynska A, Szweykowska-Kulinska Z. MicroRNA-mediated regulation of flower development in grasses. Acta Biochim Pol. 2016;63(4):687–692.
Jiang G, Xiang Y, Zhao J, Yin D, Zhao X, Zhu L, Zhai W. Regulation of inflorescence branch development in rice through a novel pathway involving the pentatricopeptide repeat protein sped1-D. Genetics. 2014;197(4):1395–1407.
Song G, Zhang R, Zhang S, Li Y, Gao J, Han X, Chen M, Wang J, Li W, Li G. Response of microRNAs to cold treatment in the young spikes of common wheat. BMC Genomics. 2017;18(1):212.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–360.
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–295.
Wang G, Yin H, Li B, Yu C, Wang F, Xu X, Cao J, Bao Y, Wang L, Abbasi AA et al. Characterization and identification of long non-coding RNAs based on feature relationship. Bioinformatics. 2019;35(17):2949–2956.
Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, Mitchell AL, Potter SC, Punta M, Qureshi M, Sangrador-Vegas A et al. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 2016;44(D1):D279-285.
Krzywinski M, Schein J, Birol I, Connors J, Gascoyne R, Horsman D, Jones SJ, Marra MA. Circos: an information aesthetic for comparative genomics. Genome Res. 2009;19(9):1639–1645.
Chen C, Chen H, Zhang Y, Thomas HR, Frank MH, He Y, Xia R. TBtools: an integrative toolkit developed for interactive analyses of big biological data. Mol Plant. 2020;13(8):1194–202.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001;25(4):402–408.
Wang Y, Gao L, Li J, Zhu B, Zhu H, Luo Y, Wang Q, Zuo J. Analysis of long-non-coding RNAs associated with ethylene in tomato. Gene. 2018;674:151–160.
Li J, Ma W, Zeng P, Wang J, Geng B, Yang J, Cui Q. LncTar: a tool for predicting the RNA targets of long noncoding RNAs. Brief Bioinform. 2015;16(5):806–812.
We thank Ms. Xiuping Xu for technical assistance on scanning electron microscopy.
This work was supported by the National Key Research and Development Program of China (2016YFD0101004), the Ministry of Agriculture of China (2016ZX08009-003), and the National Natural Science Foundation of China (31500298).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Morphology of wheat spike at six developmental stages. a-f Morphology of spikes at the six developmental stages (S1-S6) under a stereomicroscope. g-l The scanning electron microscopy images of the spikes described in a-f. SM, spikelet meristem. GP, glume primordium. LP, lemma primordium. FM, floret meristem. STP, stamen primordium. PP, pistil primordium. ANP, anther primordium. AWP, awn primordium. Bars: a-f, 1mm; g, 0.1 mm; h, 0.2 mm; i, 0.3 mm; j, 0.4 mm; k, 0.3 mm; l, 0.5 mm. Supplementary Figure S2 Venn diagrams of the lncRNAs expressed in wheat spike. Six developmental stages (S1 to S6) are indicated by different colors. Supplementary Figure S3 Numbers of DE lncRNAs among developmental stages. The numbers of DE lncRNAs with increased or decreased expression between two compared stages are shown in red or blue, respectively. DE lncRNAs were filtered according to p-value < 0.05 and log2 (fold change) > 1 or < − 1. Supplementary Figure S4 Quantitative RT-PCR analysis of six representative lncRNAs. The expression levels of six representative lncRNAs are shown based on the data from qRT-PCR (blue) and RNAseq (red), respectively. The wheat ACTIN was used as an internal reference to normalize the qRT-PCR results. Data are from three biological replicates and error bars indicate SD.
RNA-seq reads and genome mapping rates. Supplementary Table S2 Coding genes identified in wheat spike. Supplementary Table S3 Differentially expressed coding genes among six developmental stages. Supplementary Table S4 LncRNAs identified in wheat spike. Supplementary Table S5 LncRNAs expressed at individual developmental stage. Supplementary Table S6 Expression level and k-means clustering assignment of differentially expressed lncRNAs. Supplementary Table S7 DE lncRNA-mRNA pairs in cis- or trans-regulatory mode. Supplementary Table S8 Stage-specific lncRNAs and their target coding genes. Supplementary Table S9 LncRNAs as putative MiRNA precursors. Supplementary Table S10 LncRNAs as putative miRNA target mimics. Supplementary Table S11 Primers used in this study.
About this article
Cite this article
Cao, P., Fan, W., Li, P. et al. Genome-wide profiling of long noncoding RNAs involved in wheat spike development. BMC Genomics 22, 493 (2021). https://doi.org/10.1186/s12864-021-07851-4