Genome-wide discovery and characterization of maize long non-coding RNAs

1While the central dogma defines the primary role for RNA as a messenger molecule in the process of gene expression, there is ample evidence for additional functions of RNA molecules. These RNA molecules include small nuclear RNAs (snRNAs), small nucleolar RNAs (snoRNAs; mainly tRNAs and rRNAs), signal recognition particle (7SL/SRP) RNAs, microRNAs (miRNAs), small interfering RNAs (siRNAs), piwi RNAs (piRNAs) and trans-acting siRNAs (ta-siRNAs), natural cis-acting siRNAs and long noncoding RNAs (lncRNAs). lncRNAs have been arbitrarily defined as non-protein coding RNAs more than 200 bp in length, distinguishing them from short noncoding RNAs such as miRNAs and siRNAs. Rather, lncRNAs have been reported to influence the expression of other genes. Based on the anatomical properties of their gene loci, lncRNAs were further grouped into antisense lncRNAs, intronic lncRNAs, overlapping lncRNAs that in part overlap protein-coding genes and intergenic lncRNAs. lncRNAs are usually expressed at low levels, lack conservation among species and often exhibit tissue-specific/cell-specific expression patterns.

2With the advent of genomic sequencing techniques, genome-wide scans for lncRNAs have been conducted via cDNA/EST in silico mining, whole genome tilling array and RNA-seq approaches and epigenetic signature-based methods. Thousands of lncRNAs have been identified in a number of species. For example, approximately 10,000 human lncRNAs were uncovered by the ENCODE Project. The finding that several hundred human lncRNAs interact with chromatin remodeling complexes suggests that they have functional significance. Indeed, some lncRNAs have been shown to influence human disease, plant development, and other biological processes.

3Although less well characterized than mammalian lncRNAs, plant lncRNAs have defined functional roles. Vernalization in Arabidopsis is influenced by lncRNAs COOLAIR (an antisense lncRNA) and COLDAIR (an intronic lncRNA). INDUCED BY PHOSPHATE STARVATION1 is a member of the TPS1/Mt4 gene family that acts as a miR399 target mimic in fine tuning of PHO2 (encoding an E2 ubiquintin conjugase-related enzyme) expression and phosphate uptake in Arabidopsis, tomato and Medicago truncatula but does not encode a protein. Enod40 was also identified as a lncRNA involved in nodulation. Genome-wide scans for lncRNAs have also been performed in Arabidopsis thaliana, Medicago truncatula, Oryza sativa and Zea mays. In maize, an in silico bioinformatic pipeline was used on a limited set of full-length cDNA sequences to identify 1,802 lncRNAs, of which 60% are likely to be precursors of small RNAs. Each of the lncRNA surveys in plants has uncovered a substantial number of lncRNAs, which are often expressed at low levels in a tissue-specific manner as in humans and other mammals, and act as natural miRNA target mimics, chromatin modifiers or molecular cargo for protein re-localization.

4To identify a more comprehensive set of maize lncRNAs, we integrated the information from available public ESTs, maize whole genome sequence annotation, and RNA-seq datasets from 30 different experiments and developmental stages in the reference genotype of maize-B73. In total, 1,704 high-confidence lncRNAs (HC-lncRNAs) and 18,459 pre-lncRNAs (which are likely to be precursors of small RNAs) were identified in this analysis. The expression patterns and potential regulatory roles of these lncRNAs were examined in 30 B73 experiments and at several well-characterized loci. Finally, we explored the regulatory variation of lncRNAs in an RNA-seq dataset of shoot apices from 105 genotypes of the maize intermated B73 × Mo17 recombinant inbred line (IBM-RIL) population to map the genetic factors underlying the expression variation of lncRNAs. These expression quantitative trait locus (eQTL) mapping results enhance our understanding of the inheritance of lncRNA expression in plants.

5We sought to identify a relatively comprehensive set of maize lncRNAs. To achieve this, it is important to remove potential pseudogenes that have acquired nonsense or missense mutations as well as non-coding RNA precursors that will give rise to known classes of RNAs such as tRNAs, rRNA, and snRNAs. A comprehensive set of transcripts for the reference genotype B73 was developed by combining data from two sources: the maize working gene set transcripts; and de novo transcript assemblies from RNA-seq datasets from 30 different experiments (Figure 1A). There are 110,028 loci (136,774 transcript isoforms) in the working gene set (WGS) of the maize genome annotation. This set of genes consists of both computational predictions of genes as well as EST collections from a variety of tissues. Many analyses in maize utilize the 39,656 genes in the filtered gene set (FGS), a subset of the WGS that was selected based upon sequence similarity to other species and the existence of putative full-length coding sequences. However, the WGS may include lncRNAs. We also developed a set of transcript assemblies based upon 806 million uniquely mapped reads from 30 different experiments of the reference genotype-B73. These sequences were used to perform de novo transcript assembly with Cufflinks and resulted in 83,623 expressed loci with 98,444 transcript isoforms, of which 16,759 loci and 17,696 transcript isoforms are not present in the WGS. The 110,028 loci (136,774 transcript isoforms) from the WGS and 83,623 loci (98,444 transcript isoforms) from the de novo transcript assemblies were combined to generate a non-redundant set of 126,787 transcribed loci (154,470 transcript isoforms) (Figure 1B,C).

6This comprehensive set of transcribed sequences from B73 was analyzed to identify putative lncRNAs. There are 33,565 loci (38,967 transcript isoforms) that are at least 200 bp in length and do not encode an ORF of more than 100 amino acids. These sequences were filtered by comparing with the Swiss-Protein database to eliminate transcripts that contain sequence similarity (E-value ≤0.001) to known protein domains. Further filtering was performed using the Coding Potential Calculator, which assesses the quality, completeness and sequence similarity of potential ORFs to proteins in the NCBI protein database. After applying these criteria, we identified 19,608 loci (20,163 transcript isoforms; Additional file 1) that encode transcripts of >200 bp but that have little evidence for coding potential, and that were considered as putative lncRNAs. These include 12,431 loci (12,647 isoforms) from the WGS and 7,177 loci (7,515 isoforms) from the de novo transcript assemblies. This set of putative lncRNAs also includes 1,580 sequences previously identified by Boerner and McGinnis.

7These 20,163 putative lncRNAs may contain precursors to small RNA molecules, such as miRNAs, short hairpin RNAs (shRNAs) and siRNAs. The putative lncRNAs were compared to a comprehensive set of small RNAs from different tissues and small RNA related mutants. More than 90% (18,459) of the putative lncRNAs have sequence similarity with small RNAs and were classified as pre-lncRNAs (Additional file 2; Materials and methods). A set of 1,704 lncRNAs that do not have sequence similarity to known classes of noncoding RNAs were defined as HC-lncRNAs (Additional file 3). These 1,704 HC-lncRNAs include 479 sequences from the WGS and 1,225 sequences from the de novo transcript assemblies (Figure 1D). The HC-lncRNAs also contain 201 (35%) of the 572 HC-lncRNAs previously identified by Boerner and McGinnis. RT-PCR was used to validate the expression and sequence for 24 lncRNAs (Figure 2). The 24 putative lncRNAs selected for validation include 18 that were present in the working gene set from the maize genome project and 6 that are novel transcripts from our assembly of RNA-seq data. RT-PCR was performed for root, leaf and shoot tissue of 2-week old B73 seedlings and the expected products were recovered for 23 of the 24 lncRNAs tested. In some cases, there was evidence for tissue-specific expression while many of the lncRNAs were detected in all three tissues. These RT-PCR bands and specific expression were largely consistent (90/96) with the RNA-seq data. For example, lncRNA (GRMZM2G549431_T01) was not detected by both RNA-seq and RT-PCR in the leaf sample. Two of the lncRNAs (GRMZM2G010274_T01 and GRMZM2G518002_T01) showed additional isoforms in some of the tissues that may reflect tissue-specific splicing variants. RT-PCR products from 10 lncRNAs were sequenced and all 10 exhibited the appropriate sequence. We proceeded to analyze characteristics, diversity and inheritance patterns of these maize HC-lncRNAs.

8A substantial number (74%) of the pre-lncRNAs have sequence similarity to repetitive sequences of maize (Figure 1E). In contrast, the majority (98%) of the HC-lncRNAs do not contain maize repetitive sequences. Taken together, over 68% (13,811) of 20,163 maize putative lncRNAs are repetitive sequences (or transposons), which is similar to the proportion of lncRNAs in mammals. While the pericentromeric regions of most maize chromosomes have lower gene densities than chromosome ‘arms’, maize lncRNAs are more evenly distributed across chromosomes (Figure 3A). The HC-lncRNAs were characterized according to the locations relative to the nearest protein-coding genes. The majority of lncRNAs (93%) are located in intergenic regions and only 7% of the lncRNAs overlap with gene sequences. Among the intergenic HC-lncRNAs, 66 (3.9%) and 209 (12.3%) are located within 5 kb upstream and downstream of genes, respectively (Figure 3B). The remaining 83.8% of intergenic lncRNAs are at least 5 kb away from the nearest gene. This proportion (83.8%) is significantly (P = 8.1E-09) higher than the proportion of FGS genes located at least 5 kb from other FGS genes (32.6%). The majority of the lncRNAs are relatively short with very few (3%) greater than 1 kb in length (Figure 3C). Most (81%) of the lncRNAs consist of a single exon (Figure 3D). While we could not directly distinguish the transcript orientation using the non-strand-specific RNA-seq, transcript orientations could be determined using the intron splicing ‘GT-AG rule’ for those HC-lncRNA genes that contain an intron. Of the 323 lncRNAs that could be oriented based on the GT-AG intron splice sites, 23 (7%) consist of antisense transcripts.

9The lncRNA sequences were compared with genomic sequences from Arabidopsis, rice and sorghum to determine the portion of lncRNAs that had similarity (BLASTN E < 1.0E-10) to these species (Figure 3E). As expected, the conservation of lncRNAs is substantially lower than that of protein coding genes in comparisons with all three species. Permutations of random samplings of intergenic or intronic DNA were used to assess whether lncRNAs exhibit the same levels of conservation for these sequences among species. The lncRNAs have sequence similarity at the same rate as observed for intergenic sequences in all three cross-species comparisons. The maize lncRNAs exhibit the same level of conservation in Arabidopsis and rice as intronic sequences (P > 0.05) but they are significantly less conserved in sorghum (P < 0.01) than are randomly selected repeat-masked intronic sequences with similar length distribution to lncRNAs (Figure 3E).

10The level of DNA methylation within and surrounding lncRNA genes was compared with that of the FGS genes in the reference genotype B73 (Additional file 4). Similar levels of DNA methylation are observed in regions 1 kb upstream and downstream of lncRNAs and FGS genes. For both the lncRNAs and FGS genes the level of DNA methylation is reduced near the transcription start and stop sites. FGS genes show substantial levels of gene body methylation in CG and CHG contexts while the gene bodies of lncRNAs do not. Gene body methylation is often associated with genes with moderate to high levels of constitutive expression and the lack of gene body methylation for lncRNAs may reflect lower or more variable expression for these genes. The CHH DNA methylation level is quite low for both FGS and lncRNA sequences.

11The tissue-specificity of lncRNA expression was explored using the RNA-seq data from 30 different samples of B73 that represent 13 distinct tissue types. The Shannon entropy, which ranges from zero for genes expressed in a single tissue to log2(Number of tissues) for genes expressed uniformly in all tissues considered, was employed to measure the tissue-specificity of lncRNA expression. Many (54%) of the lncRNAs were only detected in one of the tissues (with at least four RNA-seq reads detected) and 10% of the lncRNAs were detected in five or more tissues (Figure 4A,B). In contrast, only 8% of FGS genes were detected in only one tissue and 74% of FGS genes were detected in five or more tissues using the same expression criteria (Figure 4A,B). Interestingly, the male reproductive tissues (immature tassel, anther, and pollen) and embryo sac had more examples of lncRNA expression than other tissues (Figure 4B). An analysis of the maximum expression level (reads per kilobase per million reads (RPKM)) for all 13 tissues provided evidence that FGS genes tend to have higher expression than lncRNAs (Figure 4C). However, 20% of the lncRNAs had an expression of >5 RPKM in at least one tissue, indicating that many of these sequences do show at least moderate expression levels in some tissues. In any one tissue, a higher proportion of FGS genes were expressed relative to HC-lncRNAs and expressed FGS genes had significantly higher expression levels than expressed HC-lncRNAs. This tissue-specific expression for many of the lncRNAs suggests that the expression of these sequences is biologically controlled rather than simply reflecting ‘transcriptional noise’.

12H3K27me3 is a facultative heterochromatin mark that is often associated with tissue-specific regulation of gene expression. The levels of H3K27me3 (trimethylation of histone H3 lysine 27) for lncRNAs were assessed in five different tissues. There are differences in the relative abundance of H3K27me3 over lncRNAs in different tissues of maize (Figure S2A in Additional file 5). The tissue with the lowest average level of H3K27me3, immature tassel, also exhibits expression for more of the lncRNAs than the other tissues, suggesting that H3K27me3 may be involved in regulating tissue-specific expression for lncRNAs. To assess the correlation between expression and H3K27me3 for the lncRNAs, H3K27me3 levels were contrasted for the lncRNAs that are expressed or silent in each of the tissues for which H3K27me3 profiles were available for analysis (Figure S2B in Additional file 5). In each tissue, genes were classified as not expressed (RPKM = 0) or expressed (RPKM >1). In general, lncRNAs that are expressed tend to have lower levels (P < 0.001) of H3K27me3 while the lncRNAs silenced in any one tissue often have elevated H3K27me3 (Figure S2B in Additional file 5). The presence of H3K27me3 at silenced lncRNAs provides evidence for targeted regulation of the expression of these lncRNAs similar to what is observed at maize genes.

13The expression levels of HC-lncRNAs in shoot apices of 105 maize IBM-RILs were compared with the expression levels in the parental lines for the 141 HC-lncRNAs that have detectable expression (at least 4 reads/RIL) in more than 40% of the RILs. The expression patterns of these 141 HC-lncRNAs were compared with those of genes in the FGS. The analysis of expression levels in shoot apices of 105 IBM RILs provides evidence for higher levels of transgressive variation in expression levels of HC-lncRNAs than in FGS genes. The difference in expression for the RILs relative to B73 or Mo17 was compared by calculating (Expparents - μprogeny)/σprogeny, which is expected to be centered at zero if the RILs generally have expression levels similar to the parents. In general, the HC-lncRNAs tend to be expressed in the RILs at levels similar to their parents but they have larger variation relative to the parents than observed for FGS genes (Figure 5A,B). This larger variation for HC-lncRNAs than FGS genes may reflect the fact that most HC-lncRNAs have quite low expression levels. However, a targeted analysis of HC-lncRNAs and FGS genes with high expression levels (RPKM >10) revealed that even highly expressed HC-lncRNAs have larger expression variation than FGS genes (Figure 5C,D). The deviation of expression levels from that of the two parents was calculated as a vector (Figure 5E) and shows evidence for higher deviation for HC-lncRNAs than for the FGS genes (P = 2.15E-20) (Figure 5F). This difference between HC-lncRNAs and FGS genes is observed for highly expressed transcripts but is not detected in transcripts with differential expression between the parents.

14The expression data from the 105 IBM RILs was used to map the regulatory regions of HC-lncRNA expression. eQTL mapping was conducted for 74 HC-lncRNAs detected in at least 80% of maize RILs using the expression levels in the 105 RILs as expression traits and a set of 7,865 high quality SNP markers. A total of 72 eQTLs (α = 0.05) with a threshold logarithm of odds (LOD) ≥4.17 were identified for 49 HC-lncRNAs. The 72 eQTLs include 21 (29%) cis-eQTLs and 51 (71%) trans-eQTLs (Figure 6A; Additional file 6), of which the proportion of trans- versus cis-eQTLs is slightly higher (P = 3.21E-03) than that observed for FGS genes. Each HC-lncRNA or FGS gene was classified according to whether a higher proportion of expression variation was explained by cis- or trans-eQTL (Figure 6B). The HC-lncRNAs were more likely to have a major trans-acting eQTL than the FGS genes (Figure 6B). Previous eQTL studies in animals and plants have revealed that many loci influenced by multiple trans-eQTL have quite high levels of expression variation in segregating offspring, presumably due to the potential for segregation of multiple eQTL with different directional effects that result in transgressive segregation. This could explain why we observe higher levels of transgressive variation for lncRNAs as they are enriched for regulation by trans-eQTL relative to FGS genes (Figure 6A). The increased contribution of trans-acting regulation to expression variation for HC-lncRNAs is consistent with the observation of higher levels of transgressive segregation for HC-lncRNA expression relative to FGS gene expression.

15We also dissected the genetic factors underlying the expression variation of 67 HC-lncRNAs, which were expressed in more than 40% but less than 80% of the RILs, as these may represent HC-lncRNAs that are expressed from one haplotype but not the other. The eQTL mapping for these 67 HC-lncRNAs identified 72 eQTLs that influenced expression of 51 of these HC-lncRNAs (Additional file 7). These HC-lncRNAs are enriched for having pre-dominant effects of cis-eQTLs (72.5%) compared with the HC-lncRNAs that are expressed in over 80% of the RILs (40.8%).

16Furthermore, 460 HC-lncRNAs were expressed (with at least 4 RNA-seq reads detected) in less than 40% of the RILs. Most (80%) of these HC-lncRNAs were expressed at very low levels (the population mean is less than 5 RPKM); while only 94 HC-lncRNAs were detected with moderate expression levels (Additional file 8). Of these moderately expressed HC-lncRNAs, only six were detected in more than 10% but less than 40% of the 105 RILs. In total, 88 out of 94 moderately expressed HC-lncRNAs were detected in only one of the 105 RILs. Taken together, these results indicate that complex regulatory mechanisms may underlie HC-lncRNA expression variation.

17There are relatively few functionally characterized lncRNAs in maize. A careful analysis of the regulation of the B1 locus in maize identified a region located more than 100 kb upstream of the coding sequence that is important for regulation and paramutation of B1 expression. There is evidence for expression of a HC-lncRNA from this region that may play a role in paramutation. Similarly, we identified a HC-lncRNA (GRMZM2G580571_T01) in the regulatory region of B1, which was previously identified as required for B’ paramutation (Figure 7A). There are several other examples of maize genes with long-distance regulatory elements. The map-based cloning of a major flowering time QTL, Vegetative to generative transition1 (Vgt1), identified a conserved non-coding region located 70 kb upstream of the ZmRap2 (GRMZM2G700665) gene, which can influence flowering time. We found a HC-lncRNA (TCONS_00089485) that is expressed from the Vgt1 regulatory region (Figure 7B). This lncRNA is detected in embryo sac and ovule tissues where ZmRap2 is not detected, while ZmRap2 is expressed in other tissues where this HC-lncRNA is not detected, suggesting the potential for antagonistic expression of this HC-lncRNA and the nearby coding sequence. The cloning of a major domestication QTL in maize identified the teosinte branched1 (tb1) gene. Further analyses provided evidence for the importance of a distant enhancer located approximately 40 kb upstream of the coding sequence that may be influenced by a transposon insertion. We also identified a pre-lncRNA (TCONS_00010027) derived from this genomic region in our study. This pre-lncRNA (TCONS_00010027) has sequence similarity with small RNAs and thus may be chopped into pieces and function as a small RNA. Because lncRNAs showed strong tissue-specific expression patterns and relatively low expression levels, and none of these three lnRNAs were detected in the tissue used for the eQTL analysis with 105 maize RILs, we could not conduct eQTL mapping for these lncRNAs. The finding that lncRNAs were detected from distant regulatory regions in all three of these examples suggests that a number of the distant regulatory regions for these maize genes may potentially involve lncRNAs. The shoot apical meristem (SAM), from which all aboveground tissues of plants are derived, is critical to plant morphology and development. While SAM initiation and development is characterized by distinct transcriptional variation, we also identified a subset of putative lncRNAs exhibiting distinct expression variation during different stages of SAM ontogeny (Additional file 9). Further research will be necessary to elucidate the functional roles of maize lncRNAs.

18The advent of high-resolution tiling arrays, the emergence of new technologies in the field of RNA-seq and large-scale chromatin immunoprecipitation experiments followed by next generation sequencing (ChIP-Seq), as well as cDNA-library sequencing and serial analysis of gene expression (SAGE), have allowed the research community to quantitatively discriminate most of the cellular transcripts. Each technical advance in examining the eukaryotic transcriptome has revealed the increasing complexity of eukaryotic genome expression. One such complexity is the existence of non-protein coding genes, including short non-protein coding genes (such as small interfering RNAs and miRNAs) and long non-protein coding genes. The short noncoding RNAs are relatively well characterized and their importance in transcriptional and posttranscriptional regulation of expression of other genes is well understood. In contrast, lncRNAs have not been as comprehensively identified or studied in many plant species.

19Our analysis generated a relatively robust list of potential lncRNAs for maize. This set of lncRNAs will likely be useful for functional genomics research or the analysis of potential functional differences among maize varieties. The lncRNAs detected in this analysis were identified from analysis of RNA-seq data from a diverse set of tissues and the current WGS annotation. In total, more than 20,000 putative lncRNAs, including 1,704 HC-lncRNAs and 18,459 pre-lncRNAs, which are likely precursors of small RNAs, were identified. We have provided GTF files as supplemental tables (Additional files 2 and 3) to enable the use and display of these lncRNAs by other researchers. Our study sheds light on the features and expression inheritance patterns of lncRNAs in maize, but also complements the reference genome annotation of maize, which might further aid the functional gene cloning and trigger more comprehensive studies on gene regulation in plants.

20Despite our use of >1 billion RNA-seq reads, it is worth noting that we only detected expression for approximately 80% of the maize FGS and approximately 50% of the lncRNAs (the other half are from WGS annotations). This may indicate that a number of additional lncRNAs with tissue- or environment-specific expression have not been detected. It is worth noting that we applied an RPKM cutoff for identifying expressed lncRNAs and that most lncRNAs were expressed at relatively low levels. While caution is required when quantifying the expression levels of genes with low RNA-seq coverage, focusing the analysis on lncRNAs with moderate expression may result in loss of lncRNAs with low expression. There are several other potential limitations to our list of lncRNAs. Most of the WGS, EST/cDNAs, and RNA-seq data were obtained after the reverse transcription with polyA primers, which selected for polyadenylated transcripts, and it is possible that some lncRNAs lack poly-adenylation. We have also employed relatively strict criteria by requiring that the putative lncRNAs lack the ability to encode peptides of more than 100 amino acids or only have a weak coding potential. However, there are examples of previously characterized lncRNAs from other species that have the potential to encode peptides >100 amino acids, such as HOTAIR with 106 amino acids, XIST with 136 amino acids and KCNQ1OT with 289 amino acids. These examples are not thought to function as proteins but would not have met our relatively strict criteria for definition as lncRNAs. Although we have identified more than 20,000 lncRNAs, it is likely that additional maize lncRNAs exist and will be discovered through analysis of additional tissues and genotypes or refinement of bioinformatics methods for characterizing lncRNAs.

21As previous studies have suggested, a substantial number of lncRNAs exist in mammals and plants, and play important functional roles in human disease, plant development, and other biological processes. In this study, we integrated available transcriptome datasets in maize to identify maize lncRNAs. More than 20,000 lncRNAs were uncovered in the maize reference genome B73, of which 1,704 were considered HC-lncRNAs. These HC-lncRNAs showed similar methylation levels as protein coding genes; however, they were more likely to exhibit tissue-specific expression patterns, which were also supported by epigenetic marks. eQTL mapping of the HC-lncRNAs showed that trans-eQTL contribute more to the expression-level variation of lncRNAs. Finally, we identified lncRNAs that were derived from regulatory regions controlling Tb1, Vgt1, and B1, which are key genes of developmental and agronomic importance in maize. We present the first comprehensive annotation of lncRNAs in maize, which opens the door for future functional genomics studies and regulatory expression research. Our findings constitute a valuable genomic resource for the identification of lncRNAs underlying plant development and agronomic traits. We also identified potential genetic mechanisms that control expression variation of lncRNAs in plant genomes.

22Transcribed sequences from the maize reference inbred line B73 were collected from the Sequence Read Archive and GenBank. Data available in the Sequence Read Archive from the maize inbred line B73 included 30 RNA-seq experiments from 13 distinct tissues (leaf, immature ear, immature tassel, seed, endosperm, embryo, embryo sac, anther, ovule, pollen, silk, and root and shoot apical meristem) encompassing a total of 1.168 billion reads with read lengths ranging from 35 to 110 nucleotides (Additional file 10). The RNA-seq data were not derived from strand-specific sequencing. Hence, it was not possible to determine transcription orientation for transcripts that do not contain introns. Maize ESTs including full-length cDNAs used by Boerner and McGinnis from a vast variety of tissues and stages were also collected from GenBank and integrated with the maize B73 genome annotation (AGP v2) (Additional file 10).

23The different sequence datasets were merged into one non-redundant set of transcript isoforms in maize, which was subjected to a series of filters to eliminate potential protein-coding transcripts (Figure 1).

24For the RNA-seq data, all sequenced reads from each experiment were aligned to the maize reference genome (AGP v2) using the spliced read aligner TopHat. Then, a method of two iterations of TopHat alignments proposed by Cabili et al. was employed to maximize the use of splice site information derived from all samples. We then re-aligned each experiment using the pooled splice sites file. The transcriptome of each experiment was assembled separately using Cufflinks. To reduce transcriptional noise, only those assembled transcript isoforms that were detected in two or more experiments were retained for further analyses. Then, we compared the assembled transcript isoforms with the maize genome annotation WGS, which represents all transcript isoforms identified by the maize genome project. The RNA-seq dataset enabled us to identify 17,696 transcript isoforms from 16,759 unknown genomic loci after filtering with the WGS. For the maize genome annotation-based transcripts, we combined maize ESTs and the WGS to eliminate transcripts from the WGS that were in silico annotated without expression evidence. The non-redundant transcripts supported by ESTs and/or RNA-seq were further filtered as follows (Figure 1).

25Putative lncRNAs were arbitrarily defined as transcripts that are ≥200 bp and have no or weak protein coding ability. We used in house perl scripts to first exclude transcripts smaller than 200 bp.

26More than 95% of protein-coding genes have ORFs of more than 100 amino acids. To remove transcripts with long ORFs, which are more likely to encode proteins, a Perl script was developed to ensure that transcripts that encode ORFs of 100 or less amino acids or incomplete ORFs were considered as lncRNA candidates.

27Transcripts were aligned to the Swiss-Protein database to eliminate transcripts with potential protein-coding ability (cutoff E-value ≤0.001).

28The Coding Potential Calculator, which is based on the detection of quality, completeness, and sequence similarity of the ORF to proteins in current protein databases, was utilized to detect putative protein encoding transcripts with default parameters. Only transcripts that did not pass the protein-coding-score test were classified as lncRNAs.

29To rule out housekeeping lncRNAs (including tRNAs, snRNAs, and snoRNAs), putative lncRNAs were aligned to housekeeping lncRNA databases. The housekeeping lncRNA databases include the tRNA database downloaded from the Genomic tRNA Database; the rRNA database from the TIGR Maize Database; and the snRNAs, snoRNAs, and signal recognition particle (7SL/SRP) collected from NONCODE. lncRNA candidates that have significant (P < 1.0E-10) alignment with housekeeping lncRNAs were not included in further analyses. Small RNAs in maize, which mainly consist of miRNAs, shRNAs and siRNAs, are generated from their precursors. The small RNA precursors are a special kind of lncRNA. To uncover this kind of lncRNA, we aligned putative lncRNAs with small RNA datasets from multiple tissues, including leaf, ear, tassel, pollen, shoot and root, and different small RNA-related mutants, mop1 and rmr2, using the same cutoff values used by Boerner and McGinnis. Here, we treated the putative lncRNAs containing homologous sequences to small RNAs as likely precursors of small RNAs; however, some of them may indeed belong to lncRNAs. Conversely, although HC-lncRNAs have no significant alignment with small RNAs, they may still be precursors of small RNAs, which could be expressed at such low levels that they could not be detected using current sequencing technology. Moreover, we also annotated lncRNAs by RepeatMasker (repetitive database version 20130422 from) with default parameters. For the classification of anatomical relationships between lncRNA loci and protein-coding genes, ZmB73 5b annotation was employed to distinguish intergenic lncRNAs from neighboring protein-coding genes. The source code for the lncRNA identification pipeline was released in the GitHub Repository.

30The above protocol used to identify lncRNAs is similar to previous studies in mammals and plants. However, we employed more stringent criteria than did Boerner and McGinnis. We used ORF ≤100 amino acids as the cutoff, whereas Boerner and McGinnis used <120 amino acids, and double filters of protein-coding potential (known protein domain filter and Protein-coding-score test).

31To validate the putative lncRNAs we identified, we conducted RT-PCR of 24 putative lncRNAs in B73 and Mo17 tissues. We grew 10 plants of B73 and Mo17 and sampled the roots, leaves and shoot apices from 14-day-old seedlings. RNA from the roots, leaves and shoots of B73 and shoots of Mo17 was isolated and used for first-strand cDNA reverse transcription by ImProm-IITM Reverse Transcription System (Promega, Madison, WI, USA). A total of 24 putative lncRNAs were randomly selected for validation and RT-PCR was conducted on the lncRNAs using routine PCR programs (Tm = 60°C) with 35 amplification cycles. To control for genomic DNA contamination in our samples, the housekeeping gene Actin was used an as experimental control. All primer information can be found in Additional file 11.

32We employed 1,000 permutations of random sequences for the significance test of sequence conservation of the FGS, lncRNAs, and intergenic and intronic fragments as follows. First, we generated intronic and intergenic annotation files based on the maize WGS annotation and all transcripts identified using the RNA-seq data in this study. Second, we randomly selected a specific number (the same to that of HC-lncRNAs) of intronic and intergenic genomic regions. Third, we adjusted the selected genomic regions according to the transcript length distribution of HC-lncRNAs. Fourth, we obtained the sequences of selected genomic regions based on the repeat-masked reference genome sequence. Fifth, we aligned the randomly selected and length-adjusted intronic and intergenic sequences against the whole genomes of Arabidopsis, rice and sorghum. Sixth, we summarized the proportion of aligned selected genomic regions with the Arabidopsis, rice and sorghum genomes with a cutoff E ≤1.0E-10. Seventh, we repeated steps 2 to 6 until the total permutation number reached up to 1,000.

33Two RNA-seq datasets were collected for the analyses of variation in lncRNA expression among tissues or among different genotypes: 1) RNA-seq data from 13 distinct tissues of the inbred line B73 from 30 experiments; and 2) RNA-seq data from 2-week old seedling shoot apices of 105 maize RILs.

34The expression levels (RPKM) of all transcripts were quantified in the two RNA-seq datasets and normalized using Cufflinks v0.9.3 based on the uniquely mapped reads of each sample. The FGS genes, most of which are conserved among species and more likely to be protein-coding genes, were used as controls for the analysis. Tissue-specific analysis measured by Shannon entropy was conducted by expression-level profiling comparison between lncRNAs and the FGS. As previously reported, bisulfite sequencing was conducted on the DNA extracted from the third seedling leaf of B73. The DNA methylation levels in CG, CHG and CHH contexts of B73 were calculated for genomic regions from which lncRNAs are transcribed and the FGS and their flanking 1 kb genomic regions. H3K27me3 levels of lncRNAs and FGS were obtained from data reported by Makarevitch et al.. For comparison of epigenetic levels, transcription start and stop sites, and upstream and downstream regions were classified based on ZmB73 5b annotations.

35In our previous study, RNA-based sequencing by Illumina Hi-Seq2000 with 103 to 110 cycles were conducted on the pooled RNA samples of 2-week-old seedling shoot apices from three replicates per genotype for 105 maize intermated B73 × Mo17 recombinant inbred lines (IBM-RILs), which were derived from the cross of the inbred lines B73 and Mo17. Uniquely mapped reads were employed to quantify the expression levels of lncRNAs and the FGS. lncRNAs and genes comprising the FGS, which were detected in the IBM-RILs, were extracted for expression inheritance pattern analysis and genetic mapping. To quantify the expression inheritance of transcripts in the RILs relative to B73 or Mo17, we used a statistic calculated by (Expparents - μprogeny)/σprogeny, where Expparents shows the expression level in the two parents, μprogeny indicates the mean value of the expression level in the progeny population and σprogeny represents the standard variation of the expression level in the progeny population for a specific gene. Any specific transcript could have two adjusted values, which measure the expression level deviation from that of the two parents (Figure 5E). The higher value the statistic is, the more deviation the transcript exhibits in the progeny compared with that of the parents. This statistic is expected to be centered at zero if the RILs generally have expression levels similar to the parents.

36A high-resolution SNP genetic map of the IBM population based upon 7,856 high quality SNP markers from RNA-seq data was used to perform eQTL mapping for lncRNAs and FGS by using composite interval mapping. To obtain a global significance of 0.05 for the eQTL mapping, a permutation threshold was computed using 1,000 randomly selected e-traits × 1,000 replicates. This threshold gave a likelihood ratio test value of 19.23, which corresponds to a LOD score of 4.17 as the significant cutoff of eQTL mapping. The confidence interval of eQTL was selected based on the range of a 1.0 LOD drop on each side from the LOD peak point. If two adjacent peaks overlap in less than 10 cM, we considered them as one eQTL.

37bp: base pair; eQTL: expression quantitative trait locus; EST: expressed sequence tag; FGS: filtered gene set; H3K27me3: trimethylation of lysine 27 of histone 3; HC-lncRNA: high confidence lncRNA; IBM-RIL: intermated B73 × Mo17 recombinant inbred line; lncRNA: long noncoding RNA; LOD: logarithm of odds; miRNA: microRNA; ORF: open reading frame; PCR: polymerase chain reaction; RIL: recombinant inbred line; RPKM: reads per kilobase per million reads; SAM: shoot apical meristem; shRNA: short hairpin RNA; siRNA: small interfering RNA; snoRNA: small nucleolar RNA; SNP: single-nucleotide polymorphism; snRNA: small nuclear RNA; WGS: working gene set.

38The authors declare that they have no competing interests.

39LL, NMS, and GJM conceived the project idea. LL, SRE, NMS, JY and GJM performed data analysis. LL, RS, KP, CTY, WW, AMC, SAG, RAC, JEF, MMSE, MJS, PSS, and MCPT performed data collection. JY, PSS, MJS, and MCPT edited the manuscript. LL, NMS and GJM wrote the manuscript. All authors read and approved the final manuscript.