Population transcriptomic analysis identifies the comprehensive lncRNAs landscape of spike in wheat (Triticum aestivum L.)

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

Download PDF

Research Article

Population transcriptomic analysis identifies the comprehensive lncRNAs landscape of spike in wheat (Triticum aestivum L.)

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

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background

Long noncoding RNAs (lncRNAs) are emerging as important regulators in eukaryotes and play crucial roles in plant growth and development as well as stress response. However, current lncRNA studies were mainly performed at the individual level and those associated with spike architecture remain elusive in wheat.

Results

In this study, RNA-seq of 186 spike samples were used to identify the lncRNA landscape in wheat spike. A total of 35,913 lncRNAs as well as 1,619 lncRNA-mRNA pairs comprised of 443 lncRNAs and 464 mRNAs were obtained. The lncRNAs displayed rather low conservation among wheat and other gramineous species. Genetic diversity analysis based on population genomics showed that obvious genetic bottleneck effect occurred on these lncRNAs during wheat domestication. Regulatory network analysis revealed that 122 lncRNA could act as ceRNA to involve in endogenous competition. Association analysis of candidate lncRNA-mRNA pairs identified 170 lncRNAs and 167 target miRNAs significantly associated with spike-related traits, including lncRNA.127690.1/TraesCS2A02G518500.1 (PMEI) and lncRNA.104854.1/TraesCS6A02G050300.1 (ATG5) were related to heading date and spike length, respectively.

Conclusions

This study provided the first pan-lncRNA data of bread wheat spike, which will contribute to the future use of lncRNA for improving the spike or yield-related traits in wheat and beyond.

Wheat

Spike

lncRNA

GWAS

Population transcriptome

It has been widely demonstrated that higher eukaryotic genomes could transcribe the huge number of non-protein coding RNAs (ncRNAs), which have indispensable regulatory functions on the epigenetic control of developmental trajectories [1, 2]. ncRNAs can be classified into small RNAs, medium-sized ncRNA and long ncRNA based on their size [3]. Long ncRNA (lncRNA) is defined as transcripts with the size of longer than 200 bp and had no discernible coding potential [4], which can be further classified as long intergenic ncRNAs (lincRNAs), intronic ncRNAs (incRNAs), sense-lncRNAs and natural antisense transcripts (NATs) transcribed from the complementary DNA strand of the associated genes [4, 5]. LncRNAs were long considered little beyond transcriptional noise; however, there is increasing evidence that lncRNAs play key roles in diverse biological processes across eukaryotes [4, 6]. Generally, lncRNA was reported to function as either cis- or trans- element through affecting neighboring loci or performing distal regulatory effect, as enhancers, decoys, guides or scaffolds for protein complexes to control gene expression [3]. With the development of next generation sequencing technologies, extensive studies have conducted in genome-wide characterization and and functional validation of lncRNAs. In plants, lncRNAs have been functionally characterized as participating in multiple biological processes such as flowering time, sexual reproduction, root organogenesis, biotic and abiotic stress responses [7, 8]. For instance, 24 lncRNAs and their potential regulatory genes were identified to be related to spike development based on RNA-seq data of six wheat spike tissue [9]. Forty-eight cis-lncRNA and 31 trans-lncRNA were identified to be associated to leaf-color of Ginkgo biloba L [10]. Despite increasing efforts, the second-generation transcriptome sequencing will have a higher false positive of lncRNA identification due to its short reading length [11]. Simultaneously, compared to human and other model animals, the biogenesis, function and evolution of lncRNA in plants, especially those underlying important agricultural traits, is not well understood.

Wheat is one of the most important staple crops worldwide, accounting for approximately 30% of the global cultivated area, and providing 20% of the world's food consumption [12]. It is estimated that wheat production needs to be increased by 1.5% annually to meet the demands of 9 billion people in 2050 [13]. Consequently, improvement of wheat yields is urgently needed. Similar to other cereal crops, such as maize, rice and barley, the yield of wheat per plant is a complex trait determined by spike number per plant, grain number per spike (GN) and one-thousand kernel weight (TKW), all of which depend largely on spike architecture and are thus considered as spike-related traits. In light of their importance, extensive studies have been performed to identify the genes and loci underlying spike architecture and spike-related traits. The Q locus, encoding an AP2-like transcription factor, is one of the well-known and widely studied domestication genes regulating spike compactness, brittle rachises and grain morphology [14]. Based on the 90 K iSelect SNP genotyping assay, 306 loci were found to be significantly associated with heading and flowering dates in 375 Chinese common wheat and these loci were also associated with spike length, peduncle length, fertile spikelet number, cold resistance, and tiller number [15]. FLOWERING LOCUS C (FLC) is reported act as a major regulator gene of flowering time in Arabidopsis [16]. Vernalization will lead to the enrichment of trimethylated histone H3 Lys on the chromatin of FLC, and lead to the stable epigenetic inhibition of FLC in Arabidopsis. Meanwhile, vernalization mediated epigenetic inhibition of FLC has been confirmed to require a long intron noncoding RNA [17], indicating the important role of lncRNA in flowering stage. In addition, one lncRNA in rice, XLOC_057324, is demonstrated to play a role in panicle development and fertility using qRT-PCR, in situ hybridization and mutant plant [18]. Two hundred and thirty-seven and 20 lncRNAs were identified to related with pollen mother cell (PMC) and embryo sac mother cell (EMC) meiosis in autotetraploid rice and were verified by CRISPR/Cas9 editing and over-expression experiment [19]. In wheat, 2,753 differentially expressed lncRNAs were analyzed in six spike development tissue containing inflorescence meristem, floral meristem, spikelet meristem, glume primordium, stamen and pistil primordia, showing the regulatory networks in wheat spike development [9]. Although many researches on plants’ lncRNA have been reported, there were few studies about lncRNA that associated with wheat yield-related traits, especially at the population transcriptome rather than the individual level.

In this study, we analyzed strand-specific RNA-seq of 186 samples obtained from wheat spike at heading and flowering stage. We constructed a comprehensive pan-lncRNA dataset containing 93 different wheat genotypes. Comparative analysis between coding RNA and lncRNA emphasizes the unique function and potential mode of action. Comprehensive evolutionary analysis revealed that lncRNA may also be selected together with functional genes in the process of wheat domestication and improvement. Based genetic variation and known QTLs, functional analysis revealed agronomic traits that could be influenced by lncRNA action. Integrating with multiple aspect datasets, our results indicated that the potential mechanism and evolution of lncRNA regulating key yield-related traits of wheat.

Genome-wide identification of lncRNAs in 186 wheat spike samples

To comprehensively obtain lncRNA related to wheat spike, we used a total of 186 RNA-seq datasets of the spike samples collected from 93 bread wheat accessions to identify the lncRNA. Totally, about 4.6 Tb clean reads were aligned to the wheat reference sequence (IWGSC Refseq v1.1) and all samples showed that the average alignment rate was more than 90%. Reference guided transcriptome assembly strategies were further applied to identify lncRNAs and then 35,913 lncRNAs were obtained (Fig. 1a, Additional file 1: Table S1). Of these, a total of 4,146 (11.5%) lncRNAs were found in the plant lncRNA database (http://greenc.sequentiabiotech.com/wiki2/Main_Page) and the average sequence identity was 92.4%.

Compared to mRNA, these wheat lncRNAs showed different basic characteristics. The average of expressed samples of lncRNAs was 119, whereas that of mRNA was 161 (Fig. 1b). The median length of lncRNAs and mRNA were 1,584 bp and 3,953 bp, respectively (Fig. 1c). Meanwhile, the number of exons of lncRNA transcripts was less than protein-coding transcripts and the majority of lncRNA transcripts (89.2%) had only one exon (Fig. 1d). However, there was no significant difference of GC content between mRNAs (50.0%) and lncRNAs (49.1%) (Fig. 1e). The median expression levels of mRNAs were higher than that of lncRNAs (Fig. 1f). Furthermore, the mRNA and lncRNA density in different chromosome and subgenomes were calculated with 1Mb interval (Fig. 2a). The results showed that the whole genome had an average of 2.71 mRNAs and 2.41 lncRNAs per Mb, respectively. Meanwhile, the number of mRNA on each chromosome was moderately correlated (R = 0.40) with the number of lncRNA through regression analysis. Moreover, most lncRNAs were identified on the B subgenome (36.1 %), followed by the A (34.9 %) and D subgenome (26.6 %) (Fig. 2b). The mRNA density of D genome was significantly higher than that of A and B subgenomes (A: 2.53/Mb; B: 2.48/Mb; D: 3.33/Mb), while the opposite trend was found on lncRNAs that D subgenome displayed the lowest lncRNA density among the three subgenomes (A: 2.49/Mb; B: 2.45/Mb; D: 2.37/Mb). Among all lncRNA, long intergenic noncoding RNA (lincRNA) (30,974) has the largest abundance (Fig. 2c). Interestingly, 36 lincRNA were identified in 15 centromeric intervals, of which 2D centromeric region had the most lncRNAs. Regression analysis showed a significant correlation between lincRNA number and centromere region size (R = 0.47, P-value = 0.0054). The mean expression level (FPKM) of these lincRNA was 1.25, which was relatively lower than other lncRNAs.

Identification of coding genes potentially targeted by lncRNAs

Correlation of expression level between lncRNAs and protein-coding genes can imply involvement in common biological processes [20]. Then, a total of 1,619 lncRNA-mRNA pairs (187 cis- and 1,432 trans-pairs) were identified in 443 lncRNA and 464 mRNA based on the physical position, expression correlation and LncTar tool (Fig. 2d, Additional file 1: Table S2). The distance between cis-lncRNAs and target mRNAs ranged from 1 bp to 4,782 bp with a mean length of 503 bp (Fig. 2e). The relationship cis-lncRNAs and target mRNAs was one to one, while one-to-many relationship was mainly identified between trans-lncRNAs and target genes (Fig. 2f). Subsequently, we performed GO and KEGG enrichment analysis of the targeting genes. Results showed that 63 GO term were significant enriched (FDR < 0.05), of which 2 and 7 terms belonged to cellular component and molecular function, respectively. The top enriched term contained pollen tube development (GO:0048868, 8.90E-16), pollen tube growth (GO:0009860, 1.50E-13), and cell tip growth (GO:0009932, 3.70E-13) (Additional file 1: Table S3, Additional file 2: Figure S1a). KEGG results showed that 11 pathways were significantly enriched, of which also contained regulation of pollen tube growth pathway (Additional file 2: Figure S1b). Enrichment analysis showed that most target coding genes were related to growth and development, especially pollen growth, suggesting the potential biological function of these lncRNAs.

Comparative analysis of conserved lncRNA sequences in homolog species

To understand the evolution of lncRNA, we performed similarity alignment and conservation analysis of the lncRNA among A. thaliana, S. bicolor, Z. mays, O. sativa, H. vulgare, Ae. tauschii, T. urartu, T. turgidum and T. aestivum. Results showed that 8,133 out of 35,193 wheat lncRNAs (23.1%) possessed homologous relationships with other species (Additional file 1: Table S4). The most homologous relationships were identified between T. aestivum and Ae. tauschii, followed by H. vulgare, T. turgidum, T. urartu, Z. mays, O. sativa, S. bicolor and A. thaliana. Most of homologous lncRNA were located in intergenic region. (Fig. 3a). GO enrichment analysis of target genes of homologous lncRNAs showed that conserved lncRNAs were involved in important biological process, such as auxin biosynthetic process (GO:0009851), protein phosphorylation (GO:0006468), and recognition of pollen (GO:0048544) (Fig. 3b).

In order to reduce the impact of different identification methods on the identification of homologous lncRNA. We further investigated whether these lncRNA sequences are highly conserved across different species genome. Wheat lncRNAs were aligned to the genomes of these eight species, and a conservation score was calculated as the product between the percentage of lncRNA sequence coverage and identity. Compared with the homologous analysis between lncRNAs and other species, more hits were obtained from the alignment of wheat lncRNAs and other species genome. Of these, the most of conserved sequences were identified between T. turgidum and T. aestivum (Additional file 1: Table S5). Based on the single copy genes among the nine species, we constructed a phylogenetic evolutionary tree (Fig. 3c). Compared with other species, the closest evolutionary relationship was also found between T. turgidum and T. aestivum. In contrast, Z. mays and O. sativa, both cereal crops, were far away from T. aestivum in the evolutionary relationship of lncRNA. These results showed that fragment introgression may be occurred between hexaploid wheat and tetraploidwheat.

Population genetic analysis of lncRNAs and target mRNAs

Based on genotyping data of 261 accessions from Triticum and Aegilops genera [21] (Additional file 1: Table S6), the genetic variations of lncRNA-mRNA pairs in A1 (wild einkorn wheat, Triticum monococcum L. ssp. aegilopoides, AA), A2 (domesticated einkorn wheat, Triticum monococcum L. ssp. monococcum, AA), A3 (urartu wheat, Triticum urartu, AA), AB1 (wild emmer wheat, Triticum turgidum L. ssp. dicoccoides, AABB), AB2 (domesticated emmer wheat, Triticum turgidum L. ssp. dicoccon, AABB), AB3 (durum wheat, Triticum turgidum L. ssp. durum, AABB), ABD1 (landrace, Triticum aestivum L. ssp. aestivum, AABBDD), ABD2 (cultivar, Triticum aestivum L. ssp. aestivum, AABBDD) and D1 (diploid, Aegilops tauschii, DD) were investigated. The fixation index (F_st) and nucleotide diversity (π) were calculated. The genetic divergence between wheat subgenome and its relatives were also detected (Fig. 4a-c). For each group, the F_st value between A2 and A3 (lncRNA: 0.82; mRNA: 0.89) tended to be largest, indicating a larger degree of divergence. The lowest F_stwas observed between ABD1 and ABD2 (lncRNA: 0.0047; mRNA: 0.0049). For each subgenome, the B subgenome had the lowest degree of divergence and there was no significant difference between A and D subgenome. The π of wild species was always larger than that of cultivated species in lncRNA and target mRNAs of different populations and subgenomes. The average π value of B subgenome was larger than that of A and D subgenome (Additional file 1: Table S7). Then, SNP density of lncRNAs and target mRNAs showed that genetic bottleneck events were found in A and D subgenome (Fig. 4d). The shifting of genetic diversity also indicated more genetic bottleneck events appeared in domesticated process (Fig. 4e). During domestication and improvement, the genetic diversity continued to shrink and the large bottleneck event occurred between A1 and A2/A3. In contrast, the B and D genomes showed the bottleneck event in domestication but genetic expansion in improvement, especially between landraces and varieties. The durum wheat and hexaploid bread wheat appeared to possess more genetic diversity after improvement, indicating the breeders would select specific genetic characteristics of wheat according to their needs in different environments.

By constructing phylogenetic evolutionary tree using genetic variation of lncRNAs and target mRNAs, we found each subgroup, expecting ABD1 and ABD2, can be almost clustered in the same branch (Additional file 2: Figure S2-S4), indicating that there were potential gene flow between ABD1 and ABD2 group and even with other subgroups. Furthermore, we identified the haplotype organization and frequency of each functional lncRNA and target mRNAs in these populations based on the resequencing data (Additional file 1: Table S8). A total 192 lncRNA and 213 mRNA were found to have the genetic variations among these populations. Among them, 3 lncRNAs and 8 target mRNAs were identified to have different haplotypes in different populations. For instance, haplotype GGGGGGTTGGGGTTCCCCTTGGAACC (76.92%) was the major haplotype in the AB3 population, but it was replaced by GGGGGGTTGGGGTTCCCCTTGGAAAA in ABD1 (93.33%) and ABD2 (100%). These results showed that these lncRNAs suffered a stronger genetic bottleneck and some loci with important candidate genes may be included in the process of domestication and improvement.

Endogenous competitive role of lncRNA in gene regulatory network

The competing endogenous RNAs (ceRNA) is also the important role of lncRNAs [22, 23]. lncRNAs can act as decoys for miRNAs to competitively inhibit their interaction with target mRNAs. Therefore, we constructed a ceRNA network of wheat spike-related lncRNA, showing the interaction relationships among lncRNA, miRNA and mRNA in wheat spike (Fig. 5a).

A total of 122 lncRNAs, 102 miRNAs and 119 mRNAs were contained in ceRNA network and 238 unique triangle pairs of lncRNA-miRNA-mRNA were identified (Additional file 1: Table S9). For instance, lncRNA.2204.2 was competitive in tae-miR319-TraesCS2A02G226000.1(TaBuB) pairs, showing the significantly negative impact between the expression level of lncRNA.2204.2 and TaBuB (Fig. 5b). The 369 bp to 389 bp of lncRNA.2204.2 cDNA is the same as the 314 to 334 bp of TaBuB, which can be cleaved by tae-miR319 (Fig. 5c-e). Moreover, tae-miR9773 involved the most of interaction relationship pairs in the ceRNA network, such as tae-miR9667b, tae-miR9780, tae-miR159, tae-miR156, tae-miR1127b and tae-miR1120b, indicating that these miRNAs may play important roles.

Functional analysis of lncRNA-mRNA pairs

In order to further explore and confirm the function of lncRNA, we focused on the 1619 lncRNA-mRNA pairs and performed the GWAS analysis with 5 agronomic traits and co-localized analysis with 124 known QTLs that associated with 13 agronomic traits. The lncRNAs or mRNAs overlapped with the known QTL or within 5 Mb upstream and downstream of GWAS signal point was considered as trait candidate gene.

Then, 54 lncRNAs and 52 target mRNAs were identified to be associated with GWAS signals, of which contained 318 lncRNA-mRNA pairs (Additional file 1: Table S10). Meanwhile, there were 130 lncRNAs and 128 mRNAs associated with known QTLs (Additional file 1: Table S11). Of these, 14 lncRNAs and 13 mRNAs were shared by GWAS signals and known QTL regions. GWAS results showed that more lncRNAs and target mRNAs were associated with spike length signals, followed by heading and flowering stage (Fig. 6). Additionally, kernel width-related QTLs had the most lncRNAs and target mRNAs. For example, lncRNA.7517.1 was identified to be related to heading stage, and its expression level was negatively correlated with the target gene. TraesCS5B02G326100.1, as the target gene of lncRNA.7517.1, is an orthologue of AtSCL5 that is a transcription factor involved in plant development [24].

Further, according to the introgression region in previous study [25], we screened the overlap regions of lncRNA and introgression. Then, 85 lncRNAs and 80 mRNAs were found in introgression interval (Additional file 1: Table S12). Meanwhile, we performed fisher's exact test and indicated a significant correlation (P-value = 2.8e-14) between the introgression intervals and lncRNA-mRNA pairs. Interestingly, the distribution of lncRNA-mRNA pairs showed that most of lncRNAs and target mRNAs were enriched in chromosome 1BS (Fig. 6). Because of the differences of wheat pedigree, we suspect that the enrichment region of 1BS may be related to wheat 1BL/1RS translocation (1B1R) line. Furthermore, a total of 55 loci (lncRNA, 33; mRNA, 32) on 1B1R region was divided into two clusters in all experimental samples by analyzing expression abundance (Additional file 2: Figure S5a). The expression level of lncRNA and mRNA between non-1B1R samples (158) and 1B1R samples (28) was also significance difference (Additional file 2: Figure S5b). GO enrichment analysis showed that chloroplast envelope (GO:0009941), cytosol (GO:0005829), maintenance of DNA methylation (GO:0010216) and chloroplast stroma (GO:0009570) were significant enriched (Additional file 2: Figure S5c). Furthermore, we performed qRT-PCR analysis for five randomly selected 1B1R lines-specific lncRNAs and five non-1B1R lines-specific lncRNAs from ten samples (five 1B1R and five non-1B1R lines). The results showed that five 1B1R lines-specific lncRNAs were highly expressed in 1B1R group but low or even not expressed in non-1B1R lines, while the expression pattern of non-1B1R lines-specific lncRNAs was opposite (Additional file 2: Figure S6).

LncRNA-mRNA related to heading and flowering stage

Gene expression is often regulated by cis-acting elements on its promotor region, such as enhancer, terminator and transcription factor (TF) binding site [26]. Six SNPs were found in gene promotor region (3) and coding region (3) of TraesCS2A02G518500.1 which was co-expressed with lncRNA.127690.1. Of these, one GWAS signal that associated with wheat heading and flowering stage (AX-110948179, C to T, Chr2A:738974737, P-value = 8.36E-05) was in the upstream of TraesCS2A02G518500.1 (Fig. 7a). We identified that the C genotype of AX-110948179 (Chr2A:742131947) form the MBSI cis-acting element, while TT genotype was the opposite (Fig. 7b). MBSI (aaaAaaC(C/G)GTTA), a TF binding region, can be targeted by MYB transcription factor gene in regulating plant growth and development, metabolism, biological and abiotic stress, hormone response and other physiological processes [27, 28]. Interestingly, genotype C corresponded to a shorter interval from heading to flowering, and a higher expression level of TraesCS2A02G518500.1, while T genotype had a longer interval heading and flowering stage and a lower expression level (Fig. 7c and 7d). The phylogenetic tree based on SNP of TraesCS2A02G518500.1 showed genetic flow might occur between tetraploid wheat and hexaploidy wheat populations (Additional file 2: Figure S7). Interestingly, TraesCS2A02G518500.1 could encode a pectin methylesterase inhibitor family protein. PMEI in rice has been reported to associated with heading stage [29]. Moreover, lncRNA.127690.1 can remotely target TraesCS2A02G518500.1 with a significant positive relationship of expression level (Fig. 7e). These results supported that lncRNA.127690.1-TraesCS2A02G518500.1 pair may play a key role in heading and flowering stage of wheat.

LncRNA-mRNA related to spike length

In addition, a GWAS signal that associated with spike length (AX-110507832, Chr6A:24725411, P-value = 2.39E-05) was identified in the downstream of lncRNA.104854.1 and its cis-regulatory mRNA (TraesCS6A02G050300.1) (Fig. 8a). Meanwhile, three SNP were found in the eighth exon region of mRNA and these three SNP changed the amino acids (Fig. 8a). Then, three main haplotypes (Hap1: CTT, Hap2: CCT, Hap3: CCC) of this mRNA were identified. The spike length of Hap1 was the longest, followed by Hap2 and Hap3 (Fig. 8b and 8c). Interestingly, lncRNA.104854.1 can compete with tae-miR5175-5p to bind to TraesCS6A02G050300.1. The expression level of lncRNA.104854.1 in Hap1 was the lowest, followed by Hap2 and Hap3 (Fig. 8d). Notably, the expression level of lncRNA.104854.1 and TraesCS6A02G050300.1 showed a significant negative relationship (Fig. 8e). Moreover, TraesCS6A02G050300 is an orthologous gene of AtATG5 encoding autophagy protein. Autophagy protein has been reported to be related to plant growth and development [30, 31]. These results indicated that TraesCS6A02G050300.1 play an important role in spike growth, and lncRNA.104854.1-TraesCS6A02G050300.1 is likely to be a key regulatory module in spike length of wheat.

Although lncRNA has been reported to involve the significant phenotype in various plant species, including defense-response and development, few studies were comprehensively performed between lncRNAs and yield-traits in wheat. Meanwhile, for the single genotype data, population analysis can use the data of multi-genotype samples to make up for genotype differences and random errors, showing a more accurate overall view [32, 33]. Here, we comprehensively identified the lncRNA at the genome-wide level consisting of 93 major Chinese varieties and landraces. Based on two coding capability prediction software and Pfam database, the coding ability of all the assembled transcripts was preliminarily predicted. To confirm the accuracy of lncRNA and its expression in the population, we retained lncRNA with the expression level and non-coding types. Finally, a total 35,913 expressed lncRNA were identified, of which 11.5% were found in plant lncRNA database. The remaining lncRNAs were not found in the database, which may be mainly due to the diversity of our samples. These results illustrated the generality of our study. Comparative analysis showed that there were differences between lncRNA and protein-coding transcripts in length, exon number, isoform number, GC content and expression level, which were basically consistent with those reported in previous studies [10, 20, 34]. Consistent to previous study [20], there were also some lincRNA within centromeric regions and the number of lincRNA was related to the length of centromeric regions, indicating the potential function of lincRNA during the reproductive development of wheat.

The wheat genome has an allopolyploid structure resulting the extensive homology across chromosomes [35]. To investigate possible sequence conservation of lncRNAs, we carried out alignment between wheat lncRNA and other eight species’ lncRNA and genome sequences. Our results showed the evolutionary pattern of lncRNA in different species is basically consistent with that of the species itself. The GO enrichment analysis of homology lncRNA illustrate the basic biological process of flowering plants, explaining the potential cause for the homology of lncRNA. In addition to lncRNA sequence homology analysis, we also performed the genetic variation analysis of lncRNA, and 261 samples of Triticum and Aegilops were selected for population genetic analysis [21]. The rapid loss of genetic diversity in lncRNAs was found in the process of domestication, and it may due to the severe bottlenecks associated with the domestication and polyploidization of Triticum. However, the genetic diversity of lncRNAs increased during improvement process on B and D subgenomes, perhaps due to hybridization with its wild relatives during modern breeding processes. Overall, lncRNAs was subjected to genetic bottleneck events similar to protein coding genes in the process of wheat evolution, indicating that it has potentially important functions.

To explore the function of lncRNAs, we further integrated the multi-omics data. Several lncRNAs can regulate mRNAs by binding or interacting with their targets and lncRNA can regulate mRNA as a ceRNA. Two lncRNAs (LOC_Os10g33940 and LOC_Os02g36880) of rice have been confirmed to be target mimics of miR160 and miR164, which has been reported to be associated with abnormal flower morphology, reduced fertility, aberrant seeds and later stages of flower development [18, 36, 37]. A total of 1,619 lncRNA-mRNA pairs and 244 ceRNA events were identified in this study. Among them, miR156 and miR159 were identified and it has been reported that miR156 can target SQUAMOSA PROMOTER BINDING PROTEIN-LIKE (SPL) genes and control the flowering pathway of Arabidopsis thaliana [38]. In addition, miR159 has also been reported to target to MYB family genes in flower development [39]. Despite the miR156/SPL and miR159/MYB module was not found in our results, the tae-miR159/PP2A was identified, which may be also related to the different development stages of the samples. Nevertheless, this results also shows the possible function of tae-miR159/PP2A module in wheat flowering development. Four and 10 lncRNAs were predicted to bind to tae-miR156 and tae-miR159, respectively, which also showed 5 shared lncRNA between tae-miR156 and tae-miR159, indicating the key function of them in heading and flowering stage.

Several introgression regions associated with significant agriculture traits have been identified in crop plants, such as rice, soybean durum wheat and bread wheat [40–42]. In this study, according the physical position of lncRNA, some lncRNAs were also located in introgression regions and QTL intervals. The 1BL/1RS line of wheat was constructed using the IRS arm from rye, which has been incorporated globally into many hexaploid wheat [43, 44]. 1BL/1RS is widely considered to have a strong relationship with wheat yield, but the days to flowering of 1BL/1RS lines is delayed [45]. In our results, most lncRNAs and target genes were located on 1BS chromosome, of which 55 loci were located in 1B1R introgression interval and mainly associated with the heading stage and spikelet number per spike. We believe that this is due to the union of lncRNA of non-1B1R and 1B1R samples in the experimental population. Because the reference genome (Chinese spring) is non-1B1R line, the 1B1R sample will show a lower expression in the translocation region than the non-1B1R sample. Further, we performed qRT-PCR of lncRNAs in 1B1R and non-1B1R samples, indicating that there are indeed differences in lncRNAs on the 1BS chromosome of translocation and non-translocation lines. Therefore, our results also determined that lncRNA on wheat 1BS was also translocated during 1R translocation (Additional file 2: Figure S5 and S6). These lncRNA may be also involved the important biological function to change the head stage and spikelet number of 1B1R lines.

Furthermore, through performing the GWAS analysis, some candidate lncRNA and target mRNAs were also identified. Autophagy protein play an important role in regulating plant growth and development as well as in stress response [46]. In this study, TraesCS6A02G050300, which was orthologous to AtATG5, was associated with spike length based on GWAS signals. The haplotype GAA has the longest spike length and higher expression of TraesCS6A02G050300, supporting the important function of this haplotype in regulating the spike length of wheat. Additionally, G genotype of SNP AX-110948179 promoted the expression of TraesCS2A02G518500, which leaded the shorter heading time of wheat in our results. These two haplotypes were identified as related to important traits and were of great significance in wheat genetic improvement. Among them, lncRNA.104854.1 and lncRNA.127690.1 could also be applied as gene resources in breeding. In the whole population, eight target gene of lncRNA have different haplotype in different population group and most haplotype changes occurred in tetraploid wheat and hexaploid wheat populations, indicating there may be more selective pressure during these periods. Among them, seven genes were annotated with Pfam database, including universal stress protein family (PF00582), zinc-finger of the FCS-type (PF04570), peroxidase (PF00141), patatin-like phospholipase (PF01734), sodium/hydrogen exchanger family (PF00999), plant invertase/pectin methylesterase inhibitor (PF04043), pectinesterase (PF01095), and Peroxidase (PF00141).

In summary, we identified 35,913 lncRNAs together with 1,619 lncRNA-mRNA pairs based on population RNA-seq dataset, which presented a comprehensive lncRNA landscape in wheat spike. This study also provided the evidence that many lncRNAs suffered a genetic bottleneck event and were located in introgression intervals during evolution. Moreover, 170 lncRNAs and 167 target mRNAs were identified to be associated with yield-related traits, including endogenous competing lncRNA. These findings provided an overall view on lncRNA for wheat genetic improvement from the transcriptional regulation level.

Plant materials and data collection

A total of 93 representative wheat lines were selected from ~5,000 accessions and used in this study [47], which were grown in the Nanyang (NY) and Luoyang (LY) trial, China in the 2018-2019 cropping season. At the heading stage, three spikes from different individuals of each accession were collected and then bulked together for subsequent experiments and RNA sequencing. All samples were collected from 9 to 11 a.m. in the field. Total RNA of each sample was isolated using the Qiagen Plant Tissue RNA Isolation kit following the instructor’s protocol. RNA-seq was performed using Illumina® HiSeq X Ten platform with paired-end 150 bp. RNA-seq data has been deposited into Genome Sequence Archive (GSA) database with the accession number of PRJCA004969 (https://ngdc.cncb.ac.cn/bioproject/browse/PRJCA004969).

Agronomic trait contained spike length (SL), grain width (GW), grain number per spike (GN), spikelet number per spike (SNS) and heading stage (HD), and these traits of each accession were collected in the NY and LY trial, China in the 2018-2019 cropping season.

The re-sequencing data of Triticum and Aegilops populations was downloaded from the Genome Variation Map published by Zhou et al. [21]. The introgression region and known QTLs region were obtained from the previous study of Chen et al. [25]. The known QTLs included grain color (GC), grain yield (GY), plant height (PH), heading date (HD), kernel number per spike (KNS), kernel weight (KW), kernel length (KL), leaf spot disease (LSD), reaction to leaf rust (LR), normalized water index (NWI), spikelet number per spike (SNS), reaction to Puccinia graminis Pers (SR), sprouting, reaction to Puccinia striiformis Westend (YR) [25].

RNA-seq analysis and lncRNAs identification

RNA-seq data was aligned to the wheat genome (IWGSC Refseq v1.1) using Hisat2 v2.2.1 [48] with default parameters. Novel and reference transcripts were assembled using StringTie v2.1.4 [49]. GTF file of each sample was merged by StringTie-merge tool and expression levels of transcripts were evaluated. Those transcript with normalized FPKM larger than 0.1 across at least 20% of the samples were considered expressed in this study. StringTie results were used to identify lncRNA. Transcripts length larger than 200 bp in length were retained for predicting coding ability using CPC and LGC with default parameters, respectively [9, 50, 51]. The Pfam database was further used to validate the protein coding ability of sequence. Then, Gffcompare was used to compare these transcripts with reference transcript. Finally, transcripts that belong to I (full contained within a reference intron), J (multi-exon with at least one junction match), O (other same strand overlap with reference exons), and U (unknown or intergenic) were considered as lncRNA [4, 9]. All expressed lncRNA were retained for subsequent analysis.

Identification of homologous lncRNAs and conserved sequences

lncRNA datasets of Arabidopsis thaliana and other cereal crops (Oryza sativa, Hordeum vulgare, Zea mays and Sorghum bicolor) were retrieved from CANTATAdb 2.0 (http://cantata.amu.edu.pl/download.php) database and IC4R information (ftp://download.big.ac.cn/bigd/GEN/IC4R-2018/1-Separated-Data/1-RNA-SeqSupported/). Homologous sequences of wheat lncRNA were identified using BLASTn with E-value cut-off 1E-05.

The wheat lncRNAs were aligned with those from two other Triticum species, one Aegilops, four other cereal crops and Arabidopsis thaliana using GMAP. The best hit of each lncRNA was extracted and that with conservation score > 0.6 [52] were used for subsequently analysis. The single copy lncRNA was used to constructed the maximum likelihood species tree using MEGA-X, with the parameter set to bootstrap of 1000 and Hasegawa-Kishino-Yano model and Gamma Distributed (G).

Genetic diversity calculation and haplotype construction

Genetic variation datasets contained 261 Triticum and Aegilops accessions (wild einkorn wheat, 31; domesticated einkorn wheat, 31; urartu wheat, 29; wild emmer wheat, 28; domesticated emmer wheat, 29; durum wheat, 13; landrace of bread wheat, 45; cultivar of bread wheat, 25; Aegilops tauschii, 30). SNPs in lncRNA and mRNA were extracted based on the chromosome location using in-house Python script. Then, nucleotide diversity (π) and fixation index (F_st) were calculated using VCFtools v0.1.16 with a sliding window of 50 Kb. The haplotype of lncRNAs and target mRNAs were constructed in nine subgroups using in-house Python scripts. Haplotypes with > 50% of the total number of subgroups were considered the major haplotype

Construction of miRNA / lncRNA-associated ceRNA networks

The target genes of lncRNA can be classified into cis- or trans-. The protein coding gene within 100 kb around lncRNA is considered to be cis-regulatory gene. In contrary, the distance between trans-regulatory gene and lncRNA is more than 100 kb. All lncRNAs without ‘N’ were included in the analysis. Meanwhile, LncTar [53] was used to calculate the free energy of lncRNA-mRNA pairs. Whether cis or trans, the expression level between lncrna and target gene also needs to be significantly correlated (spearman correlation |r| > 0.9, p < 0.05). miRNA binding sites of lncRNA and mRNA were predicted using psRNATarget [54] with default parameters. The interrelationship of lncRNA, miRNA and mRNA were identified using in-house Python script. Transcript secondary structure was predicted using RNAfold.

GWAS analysis of lncRNA-mRNA pairs

Wheat 660K genotyping datasets of all accessions from pervious study was used to perform the Genome wide associated study [47]. The physical positions of the SNPs were obtained from the Triticeae Multi-omics Center website (http://202.194.139.32/). Based on the compressed MLM model [55], GWAS was conducted using TASSEL v5.0 [56]. Then, lncRNAs within 5MB upstream and downstream of trait-associated SNPs were obtained.

GO and KEGG enrichment analysis

Gene Ontology (GO) and KEGG analysis was carried out to evaluate the potential functions and metabolic pathways for each lncRNA. The over-represented GO terms were determined using the online toolkit AgriGO v2.0 [57] with the false discovery rate (FDR) < 0.05. KOBAS v3.0 [58] was used for KEGG pathways enrichment with default parameters.

Validation by qRT-PCR

The qRT-PCR was performed for 10 randomly selected lncRNAs in five 1B1R lines and five non-1B1R lines. All reactions were performed in three separate technological replicates. RNA Easy Fast Plant Tissue Kit (Tiangen, Beijing, China) was used to extract total RNA from twenty samples and RT Master Mix Perfect Real-Time kit (Takara, Dalian, China) was used to synthesize cDNA according to the manufacturer’s instruction. qPCR reaction was performed on a QuantStudioTM 7 Flex System (Thermo Fisher Scientific, USA) using SYBR® Green Premix Pro Taq HS qPCR Kit (Accurate Biology, Hunan, China) with the following thermal cycling conditions: 95 ℃ for 30 s followed by 40 cycles of 95 ℃ for 3 s, 60 ℃ for 30 s. The expression levels were calculated using the 2^−ΔΔCT method with TaActin2 as the internal reference gene. The primers used in this study were listed in the supplementary file (Additional file 1: Table S13).

LncRNAs

Long noncoding RNAs

qRT-PCR

Quantitative real-time polymerase chain reaction

ceRNA

Competing endogenous RNAs

Transcription factor

Gene ontology

KEGG

Kyoto encyclopedia of genes and genomes

Acknowledgments

We are grateful to High-Performance Computing center of Northwest A&F University for providing computational resources in this work. We would like to acknowledge the contribution of Prof. Dejun Han in the bread wheat plant materials used in this publication.

Research involving plants

All methods were performed in accordance with the relevant Chinese national guidelines.

Author’s contributions

N.X.J. and YG designed and supervised the project. Y.G. and D.P.C. collected and generated the data, and perform analysis. Y.G. and G.Q.F performed GWAS analysis. S.T.R. and P.W.Q. contributed to qRT-PCR analysis. Y.G. prepared the draft manuscript. C.L.C, L.X.Q, and N.X.J. reviewed and revised the manuscript.

Funding

This work was funded by Tang Scholar in NWSUAF, and partially supported by the National Natural Science Foundation of China (Grant No. 31771778 and 31971885) and the Programmer of Introduction Talents of Innovative Discipline to Universities (Project 111) from the State Administration of Foreign Experts Affairs, China (Grant No. #B18042).

Availability of data and materials

The datasets supporting the conclusions of this article are included within the article and its additional files. RNA sequencing data used in this study has been deposited into Genome Sequence Archive (GSA) database with the accession number of PRJCA004969 (https://ngdc.cncb.ac.cn/bioproject/browse/PRJCA004969) and other data was provided in supporting information files.

Ethics approval and consent to participate

Experimental research and field studies on plants in this study, including the collection of plant material complied with the guidelines and legislation and NWSUAF. And the bread wheat plant materials were provided by Prof. Dejun Han and we stated that we have permission from Prof. Dejun Han to collect plant samples.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Author details

¹State Key Laboratory of Crop Stress Biology in Arid Areas, College of Agronomy and Yangling Branch of China Wheat Improvement Center, Northwest A&F University, Yangling 712100, Shaanxi, China. ²Peking University Institute of Advanced Agricultural Science, Weifang261325, Shandong, China. ³College of Biological Science and Engineering, Jiangxi Agricultural University, Nanchang 330045, Jiangxi, China.

Amaral PP, Dinger ME, Mercer TR, Mattick JS. The eukaryotic genome as an RNA machine. Science. 2008;319(5871):1787–9.
Statello L, Guo CJ, Chen LL, Huarte M. Gene regulation by long non-coding RNAs and its biological functions. Nat Rev Mol Cell Biol. 2021;22(2):96–118.
Quinn JJ, Chang HY. Unique features of long non-coding RNA biogenesis and function. Nat Rev Genet. 2016;17(1):47–62.
Chekanova JA. Long non-coding RNAs and their functions in plants. Curr Opin Plant Biol. 2015;27:207–16.
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.
Wu Z, Liu X, Liu L, Deng H, Zhang J, Xu Q, Cen B, Ji A. Regulation of lncRNA expression. Cell Mol Biol Lett, 2014;19(4):561–575.
Yu Y, Zhang Y, Chen X, Chen Y. Plant Noncoding RNAs: Hidden Players in Development and Stress Responses. Annu Rev Cell Dev Biol. 2019;35:407–431.
Jha UC, Nayyar H, Jha R, Khurshid M, Zhou M, Mantri N, Siddique KHM. Long non-coding RNAs: emerging players regulating plant abiotic stress response and adaptation. BMC Plant Biol. 2020;20(1):466.
Cao P, Fan W, Li P, Hu Y. Genome-wide profiling of long noncoding RNAs involved in wheat spike development. BMC Genomics. 2021;22(1):493.
Wu Y, Guo J, Wang T, Cao F, Wang G. Transcriptional profiling of long noncoding RNAs associated with leaf-color mutation in Ginkgo biloba L. BMC Plant Biol. 2019;19(1):527.
Stark R, Grzelak M, Hadfield J. RNA sequencing: the teenage years. Nat Rev Genet. 2019;20(11):631–656.
Miransari M, Smith D. Sustainable wheat (Triticum aestivum L.) production in saline fields: a review. Crit Rev Biotechnol. 2019;39(8):999–1014.
Zhao Q, Zhang H, Wang T, Chen SX, Dai SJ. Proteomics-based investigation of salt-responsive mechanisms in plant roots. J Proteome. 2013;82:230–253.
Zhang Z, Belcram H, Gornicki P, Charles M, Just J, Huneau C, Magdelenat G, Couloux A, Samain S, Gill BS, et al. Duplication and partitioning in evolution and function of homoeologous Q loci governing domestication characters in polyploid wheat. Proc Natl Acad Sci U S A. 2011;108(46):18737–42.
Zhang X, Chen J, Yan Y, Yan X, Shi C, Zhao L, Chen F. Genome-wide association study of heading and flowering dates and construction of its prediction equation in Chinese common wheat. Theor Appl Genet. 2018;131(11):2271–2285.
Michaels SD, Amasino RM. FLOWERING LOCUS C encodes a novel MADS domain protein that acts as a repressor of flowering. Plant Cell. 1999;11(5):949–56.
Heo JB, Sung S. Vernalization-mediated epigenetic silencing by a long intronic noncoding RNA. Science. 2011;331(6013):76–9.
Zhang YC, Liao JY, Li ZY, Yu Y, Zhang JP, Li QF, Qu LH, Shu WS, Chen YQ. Genome-wide screening and functional analysis identify a large number of long noncoding RNAs involved in the sexual reproduction of rice. Genome Biol. 2014;15(12):512.
Li X, Shahid MQ, Wen M, Chen S, Yu H, Jiao Y, Lu Z, Li Y, Liu X. Global identification and analysis revealed differentially expressed lncRNAs associated with meiosis and low fertility in autotetraploid rice. BMC Plant Biol. 2020 Feb 19;20(1):82.
Zhou Y, Zhao X, Li Y, Xu J, Bi A, Kang L, Xu D, Chen H, Wang Y, Wang YG, et al. Triticum population sequencing provides insights into wheat adaptation. Nat Genet. 2020;52(12):1412–1422.
Golicz AA, Singh MB, Bhalla PL. The Long Intergenic Noncoding RNA (LincRNA) Landscape of the Soybean Genome. Plant Physiol. 2018;176(3):2133–2147.
Liu Y, Xue M, Du S, Feng W, Zhang K, Zhang L, Liu H, Jia G, Wu L, Hu X, et al. Competitive endogenous RNA is an intrinsic component of EMT regulatory circuits and modulates EMT. Nat Commun. 2019;10(1):1637.
Zhang P, Wu W, Chen Q, Chen M. Non-Coding RNAs and their Integrated Networks. J Integr Bioinform. 2019;16(3):null.
Lee MH, Kim B, Song SK, Heo JO, Yu NI, Lee SA, Kim M, Kim DG, Sohn SO, Lim CE, et al. Large-scale analysis of the GRAS gene family in Arabidopsis thaliana. Plant Mol Biol. 2008;67(6):659–70.
Cheng H, Liu J, Wen J, Nie X, Xu L, Chen N, Li Z, Wang Q, Zheng Z, Li M, et al. Frequent intra- and inter-species introgression shapes the landscape of genetic variation in bread wheat. Genome Biol. 2019;20(1):136.
Inukai S, Kock KH, Bulyk ML. Transcription factor-DNA binding: beyond binding site motifs. Curr Opin Genet Dev. 2017;43:110–119.
Solano R, Nieto C, Avila J, Cañas L, Diaz I, Paz-Ares J. Dual DNA binding specificity of a petal epidermis-specific MYB transcription factor (MYB.Ph3) from Petunia hybrida. EMBO J. 1995;14(8):1773–84.
Dubos C, Stracke R, Grotewold E, Weisshaar B, Martin C, Lepiniec L. MYB transcription factors in Arabidopsis. Trends Plant Sci. 2010;15(10):573–81.
Nguyen HP, Jeong HY, Jeon SH, Kim D, Lee C. Rice pectin methylesterase inhibitor28 (OsPMEI28) encodes a functional PMEI and its overexpression results in a dwarf phenotype through increased pectin methylesterification levels. J Plant Physiol. 2017;208:17–25.
Thompson AR, Doelling JH, Suttangkakul A, Vierstra RD. Autophagic nutrient recycling in Arabidopsis directed by the ATG8 and ATG12 conjugation pathways. Plant Physiol. 2005;138(4):2097–110.
Yoshimoto K, Shibata M, Kondo M, Oikawa K, Sato M, Toyooka K, Shirasu K, Nishimura M, Ohsumi Y. Organ-specific quality control of plant peroxisomes is mediated by autophagy. J Cell Sci. 2014;127(Pt 6):1161–8.
Gusev A, Ko A, Shi H, Bhatia G, Chung W, Penninx BW, Jansen R, de Geus EJ, Boomsma DI, Wright FA, et al. Integrative approaches for large-scale transcriptome-wide association studies. Nat Genet. 2016;48(3):245–52.
Wainberg M, Sinnott-Armstrong N, Mancuso N, Barbeira AN, Knowles DA, Golan D, Ermel R, Ruusalepp A, Quertermous T, Hao K, et al. Opportunities and challenges for transcriptome-wide association studies. Nat Genet. 2019;51(4):592–599.
Xue L, Sun M, Wu Z, Yu L, Yu Q, Tang Y, Jiang F. LncRNA regulates tomato fruit cracking by coordinating gene expression via a hormone-redox-cell wall network. BMC Plant Biol. 2020;20(1):162.
Feldman M, Levy AA. Genome evolution due to allopolyploidization in wheat. Genetics. 2012;192(3):763–74.
Baker CC, Sieber P, Wellmer F, Meyerowitz EM. The early extra petals1 mutant uncovers a role for microRNA miR164c in regulating petal number in Arabidopsis. Curr Biol. 2005;15(4):303–15.
Liu X, Huang J, Wang Y, Khanna K, Xie Z, Owen HA, Zhao D. The role of floral organs in carpels, an Arabidopsis loss-of-function mutation in MicroRNA160a, in organogenesis and the mechanism regulating its expression. Plant J. 2010;62(3):416–28.
Yu N, Cai WJ, Wang S, Shan CM, Wang LJ, Chen XY. Temporal control of trichome distribution by microRNA156-targeted SPL genes in Arabidopsis thaliana. Plant Cell. 2010;22(7):2322–35.
Guo C, Xu Y, Shi M, Lai Y, Wu X, Wang H, Zhu Z, Poethig RS, Wu G. Repression of miR156 by miR159 Regulates the Timing of the Juvenile-to-Adult Transition in Arabidopsis. Plant Cell. 2017;29(6):1293–1304.
Han C, Zhang P, Ryan PR, Rathjen TM, Yan Z, Delhaize E. Introgression of genes from bread wheat enhances the aluminium tolerance of durum wheat. Theor Appl Genet. 2016;129(4):729–739.
Wang X, Chen L, Ma J. Genomic introgression through interspecific hybridization counteracts genetic bottleneck during soybean domestication. Genome Biol. 2019;20(1):22.
Chen Z, Li X, Lu H, Gao Q, Du H, Peng H, Qin P, Liang C. Genomic atlases of introgression and differentiation reveal breeding footprints in Chinese cultivated rice. J Genet Genomics. 2020;47(10):637–649.
Lelley T, Eder C, Grausgruber H. Influence of 1BL.1RS wheat-rye chromosome translocation on genotype by environment interaction. J Cereal Sci. 2004;39(3):313–320.
Ren T-H, Yang Z-J, Yan B-J, Zhang H-Q, Fu S-L, Ren Z-L. (2009) Development and characterization of a new 1BL.1RS translocation line with resistance to stripe rust and powdery mildew of wheat. Euphytica, 2009;169:207–213.
Qi W, Tang Y, Zhu W, Li D, Diao C, Xu L, Zeng J, Wang Y, Fan X, Sha L, et al. Molecular cytogenetic characterization of a new wheat-rye 1BL•1RS translocation line expressing superior stripe rust resistance and enhanced grain yield. Planta. 2016;244(2):405–16.
Su T, Li X, Yang M, Shao Q, Zhao Y, Ma C, Wang P. Autophagy: An Intracellular Degradation Pathway Regulating Plant Survival and Stress Response. Front Plant Sci. 2020;11:164.
Wu J, Yu R, Wang H, Zhou C, Huang S, Jiao H, Yu S, Nie X, Wang Q, Liu S, et al. A large-scale genomic association analysis identifies the candidate causal genes conferring stripe rust resistance under multiple field environments. Plant Biotechnol J. 2021;19(1):177–191.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.
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.
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.
Zhao T, Tao X, Feng S, Wang L, Hong H, Ma W, Shang G, Guo S, He Y, Zhou B, et al. LncRNAs in polyploid cotton interspecific hybrids are derived from transposon neofunctionalization. Genome Biol. 2018;19(1):195.
Lin X, Lin W, Ku YS, Wong FL, Li MW, Lam HM, Ngai SM, Chan TF. Analysis of Soybean Long Non-Coding RNAs Reveals a Subset of Small Peptide-Coding Transcripts. Plant Physiol. 2020;182(3):1359–1374.
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–12.
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 Z, Ersoz E, Lai CQ, Todhunter RJ, Tiwari HK, Gore MA, Bradbury PJ, Yu J, Arnett DK, Ordovas JM, Buckler ES. Mixed linear model approach adapted for genome-wide association studies. Nat Genet. 2010;42(4):355–60.
Bradbury PJ, Zhang Z, Kroon DE, Casstevens TM, Ramdoss Y, Buckler ES. TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics. 2007;23(19):2633–5.
Tian T, Liu Y, Yan H, You Q, Yi X, Du Z, Xu W, Su Z. agriGO v2.0: a GO analysis toolkit for the agricultural community, 2017 update. Nucleic Acids Res. 2017;45(W1):W122-W129.
Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, Kong L, Gao G, Li CY, Wei L. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011;39(Web Server issue):W316-22.

No competing interests reported.

Download PDF

Editorial decision: Major revision
01 Aug, 2022
Reviews received at journal
31 Jul, 2022
Reviewers agreed at journal
24 Jul, 2022
Reviews received at journal
21 Jul, 2022
Reviewers agreed at journal
28 Jun, 2022
Reviews received at journal
15 Jun, 2022
Reviewers agreed at journal
13 Jun, 2022
Reviewers invited by journal
13 Jun, 2022
Editor assigned by journal
13 Jun, 2022
Editor invited by journal
13 Jun, 2022
Submission checks completed at journal
13 Jun, 2022
First submitted to journal
07 Jun, 2022

You are reading this latest preprint version

Population transcriptomic analysis identifies the comprehensive lncRNAs landscape of spike in wheat (Triticum aestivum L.)

Status:

Version 1

Abstract

Background

Results

Conclusions

Figures

Background

Results

Discussion

Conclusions

Methods

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1