Systematic Identification of Methyl Jasmonate-Responsive Long Noncoding RNAs and Their Nearby Coding Genes Unveils Their Potential Defence Roles in Tobacco BY-2 Cells

Long noncoding RNAs (lncRNAs) are distributed in various species and play critical roles in plant growth, development, and defence against stimuli. However, the lncRNA response to methyl jasmonate (MeJA) treatment has not been well characterized in Nicotiana tabacum Bright Yellow-2 (BY-2) cells, and their roles in plant defence remain elusive. Here, 7848 reliably expressed lncRNAs were identified in BY-2 cells, of which 629 differentially expressed (DE) lncRNAs were characterized as MeJA-responsive lncRNAs. The lncRNAs in BY-2 cells had a strong genus specificity in Nicotiana. The combined analysis of the cis-regulated lncRNAs and their target genes revealed the potential up- and downregulated target genes that are responsible for different biological functions and metabolic patterns. In addition, some lncRNAs for response-associated target genes might be involved in plant defence and stress resistance via their MeJA- and defence-related cis-regulatory elements. Moreover, some MeJA-responsive lncRNA target genes were related to quinolinate phosphoribosyltransferase, lipoxygenases, and endopeptidase inhibitors, which may contribute to nicotine synthesis and disease and insect resistance, indicating that MeJA-responsive lncRNAs regulate nicotine biosynthesis and disease resistance by regulating their potential target genes in BY-2 cells. Therefore, our results provide more targets for genetically engineering the nicotine content and plant defence in tobacco plants.


Introduction
With the development of tiling arrays and next-generation sequencing (NGS), noncoding RNAs (ncRNAs) are being recognized as an increasingly important part of the transcriptomes of eukaryotic organisms [1,2]. Small RNAs [18-30 nucleotides (nt)], medium-sized ncRNAs (31-200 nt), and long ncRNAs (>200 nt) are the three major types of ncRNAs [3]. Among them, long noncoding RNAs (lncRNAs) have recently emerged as powerful and essential function regulators, notably in mammals [4][5][6]. However, research on diverse and complex lncRNAs in plants is still in its early stages. Small interfering RNA (siRNA) precursors, scaffolds, and molecular sponges are several functions of plant lncRNAs induced by various environmental conditions [3,7]. Some plant lncRNAs, such as COLDAIR and COOLAIR, can modulate transcription via chromatin changes and become scaffolds of DNA or protein [8,9]. Several lncRNAs function as decoys that alter target protein behaviour. During developmental transitions, Arabidopsis ASCO-lncRNA competes with NUCLEAR SPECKLE RNA-BINDING PROTEINS (NSRs) to modulate alternative splicing and gene After the read alignment, transcript assembly, and systematic exploration of the lncRNAs, we identified 7848 lncRNAs derived from 7560 loci based on RNA-seq datasets from MeJA-untreated and MeJA-treated BY-2 cells ( Figure 1A and Supplementary Table S1). It should be noted that each group of the reading segment data matched the tobacco genome at a rate of approximately 95%, and all of the lncRNAs that were eventually examined had realistic FPKM values. We know that the length of lncRNAs is generally greater than 200 nucleotides (nt), and therefore, we focused on the length distribution of the lncRNAs. In the BY-2 cells, the lncRNA transcript lengths varied from 201 to 5132 nt, with a mean of 631 nt. We found that lncRNAs of 400-599 nt were the most abundant (n = 2510), while the number of lncRNAs ranging from 200-399 nt and 600-799 nt was quite close ( Figure 1B and Supplementary Table S2). Regarding the exon count, 6971 lncRNAs had one exon, and 975 lncRNAs had two exons. The presence of multiple exons (from three to eight) was not the primary feature of tobacco BY-2 cell lncRNAs ( Figure 1C and Supplementary Table S3). Furthermore, just one lncRNA contained eight exons, the most exons of any lncRNA. Then, we divided the lncRNAs into four primary types based on their location in the genome with neighbouring genes ( Figure 1D). In detail, there were 7832 intergenic transcripts, nine transcripts of generic exonic overlap with reference transcripts, five transcripts of exonic overlap with reference on the opposite strand, and two transcripts of transfrag falling entirely within the reference intron. Therefore, the intergenic lncRNAs appeared to make up the majority of the lncRNAs (99.796%), which had an average GC content of 0.39, close to the genomic sequence of N. tabacum (0.39) ( Figure 1D).

The lncRNAs in BY-2 Cells Have Extreme Genus Specificity for Nicotiana
To comprehensively investigate the evolutionary relationship of lncRNAs in tobacco BY-2 cells, we retrieved the lncRNAs and phylogenetic relationships of thirty-nine species, including eudicots (twenty), monocots (fourteen), basal angiosperms (one), ferns (one), mosses (one), and green algae (one) ( Figure 2B). Using tobacco BY-2 cells as an example, only 4580 (58.359%) of the 7848 lncRNAs were homologous or paralogous with those of the above species (Figure 2A and Supplementary Tables S4 and S5). Within the four Solanaceae species, there was a relatively high sequence similarity, with 6535 (83.269%), 3921 (49.962%), 9 (0.115%), and 10 (0.127%) lncRNAs in N. tabacum, N. benthamiana, S. tuberosum, and S. lycopersicum having homology with tobacco BY-2 cells, respectively. It should be noted that the BY-2 cells were derived from N. tabacum, and therefore 6535 lncRNAs were paralogous or the same lncRNAs. Furthermore, we found that 822 lncRNAs of tobacco BY-2 cells were homologous between N. tabacum and N. benthamiana, indicating that they might be highly conserved in Nicotiana ( Figure 2C). The lncRNAs in the BY-2 cells were approximately one-fourth of all the lncRNAs in N. tabacum (31,028). The maximum percentage of the homologous was only 0.025% in other non-Solanaceae, and no homologous lncRNAs were detected in most of the species. These findings suggested that lncRNAs in tobacco BY-2 cells were poorly conserved with most plants but had extreme genus specificity in Nicotiana.   Table S6). In addition, the expression of the lncRNAs in the MeJA-treated group varied from 0 to 4957.32. The top three most highly expressed lncRNAs were TCONS 00015430 (FPKM = 4957.32), TCONS 00091915 (FPKM = 2790.99), and TCONS 00091916 (FPKM = 1161.31). The three lncRNAs described above were always the top three most abundant lncRNAs. According to the results of the differential expression analysis, 629 significantly differentially expressed (DE) lncRNAs were identified as MeJA-responsive lncRNAs, of which 367 were upregulated and 262 were downregulated (p value ≤ 0.01; |log 2 FC| ≥ 1) ( Figure 3B and Supplementary Table S7). Figure 3C,D depict the expression pattern heatmaps of MeJA-downregulated and -upregulated lncRNAs, respectively. Notably, all of the MeJA-responsive lncRNAs were either unknown or intergenic lncRNAs. The MeJA-upregulated lncRNAs had an upregulation fold change ranging from 5.94489 to 1.00241 (log 2 FC), with TCONS 0010583 showing the highest upregulation fold change. The downregulation fold change of the MeJA-downregulated lncRNAs varied from 3.31928 to 1.00052 (log 2 FC), with TCONS 00005455 showing the highest downregulation fold change. The expression pattern, member number, and fold change of the MeJA-responsive lncRNAs implied that upregulated lncRNAs constituted the major response group following MeJA treatment in BY-2 cells.

Potential Up-and Downregulated Target Genes of MeJA-Upregulated lncRNAs Are Responsible for Different Biological Functions and Metabolic Patterns
Because lncRNAs play vital roles in regulating gene expression, and the potential function of lncRNAs has not been previously reported in tobacco BY-2 cells, identifying and analysing their target genes might aid in our understanding of their activities. We found 9571 lncRNA-mRNA pairs, consisting of 4225 lncRNAs and 7085 adjacent coding genes within 100 kilobases of the lncRNA. Among them, 242 significantly upregulated MeJAresponsive lncRNAs and 136 downregulated MeJA-responsive lncRNAs corresponded to 518 mRNAs and 315 mRNAs via "pairs" of cis genes, respectively ( Figure 4A,B and Supplementary Table S8). In addition, three mRNAs were associated with one MeJA-upand downregulated lncRNA ( Figure 4B). Among these putative target genes of MeJAupregulated lncRNAs, 102 genes showed differential expression between the MeJA-treated and control samples, with 83 upregulated and 19 downregulated ( Figure 4C). Interestingly, these two groups of target genes were enriched in significantly different biological functions and metabolic patterns. In the biological process category, 34 GO categories were enriched among the upregulated target genes, whereas 24 GO terms were enriched among the downregulated target genes ( Figure 4D,F). In addition, metabolic and biosynthetic process genes were highly prominent, prompting us to focus on metabolic processes. Figure 4E demonstrates that the upregulated target genes for the upregulated MeJA-responsive lncRNAs were significantly enriched in the "metabolism" category (A 09100) ( Figure 4E). Moreover, the "environmental adaptation" (B 09159), "metabolism of terpenoids and polyketides" (B 09109), "metabolism of cofactors and vitamins" (B 09108), "metabolism of other amino acids" (B 09106), "amino acid metabolism" (B 09105), and "lipid metabolism" (B 09103) suggested that the upregulated target genes of the MeJAupregulated lncRNAs initiated plant defence through a string of metabolic processes. However, the downregulated genes were not strongly linked to metabolic activity, and these genes seem to play a role in genetic information processing ( Figure 4D,G). This implied that the expression levels of the genes involved in the growth and developmental processes were downregulated by the lncRNAs when the plants were stimulated.

Response-Associated Target Genes of Some Conserved MeJA-Responsive lncRNAs Are Involved in Plant Defence and Stress Resistance
The enrichments in the "response to chemical", "response to stress", "response to endogenous stimulus", "response to abiotic stimulus", "response to biotic stimulus", and "response to external stimulus" suggested that the up-and downregulated target genes of the MeJA-upregulated lncRNAs were related to diverse stimuli and plant defence. Furthermore, the response-associated up-or downregulated target genes did not share any members, and several genes worked on different response processes ( Figure 5A). We also depicted the expression patterns of all the response-associated target genes and their respective lncRNAs ( Figure 5B, Supplementary Table S9). We found a series of target genes that were significantly upregulated and associated with plant resistance to insects and diseases, such as linoleate 9S-lipoxygenase, which plays a role in plant-pathogen interactions (LOC107770253), and 21 kDa seed protein-like with endopeptidase inhibitor activity (LOC107782545) ( Figure 5) [41,42]. In Solanaceae, we found the lncRNAs homologous or paralogous to N. tabacum and N. benthamiana in the BY-2 cells and speculated that the target genes of these conserved lncRNAs might have their own specificity. Among them, the upregulated MeJA-responsive lncRNAs with upregulated target genes were TCONS_00042907, TCONS_00042744, TCONS_00037532, and TCONS_00076889, and the upregulated lncRNA with downregulated target genes was TCONS_00136591 ( Figure 5C). Specifically, TCONS_00042907 targeted LOC107829122, which encodes a quinolinate phosphoribosyltransferase (QPT) family member [43]. We discovered that the transcript levels of LOC107829122 rose approximately 50-fold and that the expression of TCONS 00042907 was detectable exclusively in the MeJA-treated group ( Figure 5B, Supplementary Table S9). In addition, LOC107826152 was the target gene of TCONS_00037532 and encoded the E3 ubiquitin-protein ligase RMA1H1-like, which has been associated with drought stress [44]. However, the high-affinity nitrate transporter 2.4-like (LOC107808057) was the only downregulated conserved target gene corresponding to the upregulated lncRNAs. In summary, the response-associated target genes of conserved MeJA-responsive lncRNAs might be involved in plant defence and stress resistance.   Table S10). As shown in Figure 6, all the MeJA-responsive lncRNAs with target genes possessed many motif elements related to the light response. The 242 downregulated MeJA-responsive lncR-NAs with target genes shared 19 similar functional cis-acting elements, with 136 downregulated MeJA-responsive lncRNAs having target genes, among which "MeJA response", "abscisic acid response", "auxin response", and "anaerobic induction" stood out ( Figure 6B,C).
It was not surprising that MeJA-associated cis-acting elements, the TGACG-motif and CGTCA-motif, appeared in the upstream sequences of the MeJA-responsive lncRNAs. Only the upregulated MeJA-responsive lncRNAs with the target genes had circadian control elements (the CAAAGATATC-motif), which was likely because the MeJA treatment caused some lncRNAs to be upregulated and interfered with the normal rhythmic activity of the plant. After analysing the response-related upregulated target genes for the upregulated MeJA-responsive lncRNAs with cis-acting elements, we found many MeJA-responsive elements upstream of TCONS_00042744 and TCONS_00042907 ( Figure 6D). The target gene of TCONS_00037532 is the E3 ubiquitin-protein ligase RMA1H1-like that is associated with plant drought stress, and the target gene of TCONS_00042744 (LOC107829039) encodes an isoform of the lipoxygenase-6 protein, which plays an essential role in the plant's response to pathogen and pest attack. Interestingly, the cis-elements of the "circadian control" and "defence and stress response" appeared upstream of TCONS_00042744, which probably conferred the disease resistance to the corresponding target genes. We hypothesized that after the MeJA treatment, the lncRNAs might be activated by their upstream MeJA-responsive cis-acting elements. Again, the presence of plant defence-related elements allowed the lncRNAs and their target genes to be involved in plant disease resistance.

Discussion
A large number of lncRNAs are associated with plant defence, but there have been no reports about them in tobacco BY-2 cells after MeJA treatment. Tobacco BY-2 cells, which were established from a callus, have been considered to be the "HeLa" cell line in higher plant cell biology [45]. MeJA has been used in experiments to mimic herbivory and induce plant defences, and therefore BY-2 cells treated with MeJA helped us to accurately study the lncRNAs involved in plant defence [28]. In this study, we systematically found a total of 7848 potential lncRNAs in tobacco BY-2 cells, using the available tobacco genome from a transcriptomic perspective (Figure 1). The number of identified lncRNAs was smaller than that identified in axillary tobacco buds (13694) and was more significant than that found in tobacco roots (7423), which might be partly due to the highly tissue-specific and cell type-specific lncRNAs [46,47]. Many studies have concluded that lncRNAs vary widely across species, tissues, and cell types, indicating that most are not conserved at the sequence level [48,49]. The low sequence identity among different plants is consistent with previous findings (Figure 2) [50]. In addition, we also found that lncRNAs in Nicotiana likely have the same evolutionary origin, suggesting that homologous lncRNAs likely share similar biological functions.
When animals or insects attack plants, the plants produce MeJA and initiate a number of plant defence mechanisms, such as deterring pest-landing, preventing feeding, and reducing palatability [15][16][17]. The MeJA-responsive lncRNAs that were significantly upregulated after the MeJA treatment were the most dominant class, which was likely regulated by MeJA-responsive cis-acting elements upstream of the lncRNAs. The TGACG motif and CGTCA motif were the main classes of the MeJA-responsive cis-elements (Supplementary Table S10). The tandem TGACG motifs can drive the expression of the 35S promoter in the roots and show responsiveness to signals, such as jasmonic acid, the auxin-related herbicide 2,4-dichlorophenoxyacetic acid, and hydrogen peroxide associated with biotic and abiotic stress [51,52]. In addition, the TGA TFs, the basic leucine zipper subfamily members associated with disease resistance and development, have been discovered to bind to the TGACG motif, and some TGA TFs from A. thaliana were shown to be necessary for the activation of jasmonic acid/ethylene-induced defence responses [53][54][55].
After ancestral hybridization, N. tabacum eventually became an allotetraploid. The biological reasons for its stable existence must be heterozygous dominance and its derivation of a more powerful gene function and regulatory network [38,39]. Previously, we have explored a series of transcriptional regulators and the JA zinc finger expressed in inflorescence meristem domain proteins associated with alkaloid formation [56,57]. It has been confirmed that lncRNAs are typically found near the genes they regulate, which is known as a cis-acting mechanism [48]. It is important to note that the lncRNAs can regulate some cis-coding genes, but they could also be coding genes that act in trans. Trans-acting lncRNAs, like mRNAs, are transcribed, processed, and then depart their transcription locations to carry out their role somewhere else [58,59]. Here, we focus only on the molecular functions that may result from cis-action. In our research, 242 DE upregulated MeJA-responsive lncRNAs, and 136 DE downregulated MeJA-responsive lncRNAs corresponding to 518 PCGs and 315 PCGs, respectively, within a genomic region of 100 kb, were coexpressed ( Figure 4, Supplementary Table S8). Furthermore, 20 of them were associated with the GO term "response to chemical, stress, endogenous stimulus, abiotic stimulus, biotic stimulus, and external stimulus" (Figure 5), suggesting that these lncRNAs may play vital roles in plant defence by modulating their target genes. Indeed, the presumably largest group of cis-acting lncRNAs are those that function to strengthen target gene expression, similar to how enhancers work [60]. The potential activation ability of most MeJA-responsive lncR-NAs on the target genes was also shown by comparing expression patterns ( Figure 5). For example, TCONS_00042907 targeted LOC107829122, which encodes a QPT family member, and the transcript levels of LOC107829122 rose approximately 50-fold. The expression of TCONS_00042907 was detectable exclusively in the MeJA-treated group ( Figure 5B, Supplementary Table S9). In topping-treated tobacco, the enhancement of nicotine biosynthesis and the upregulation of QPT2 are achieved by the trans-acting lncRNA (nta-eTMX27)mediated inhibition of the miRNA expression [61]. This result provides another lncRNA case for the regulation of QPT from the cis-acting perspective.
QPTs are the rate-limiting steps of plant nicotine synthesis, and the importance of nicotine for plant defence cannot be overstated [14,62]. In addition, we also described other target genes of MeJA-responsive lncRNAs that are relevant to plant resistance to stimuli.
The action of lipoxygenases (9-LOX and 13-LOX) or α-dioxygenase (α-DOX) initiates the upstream biosynthesis of plant JA, and LOX-derived oxylipins can be involved in the pathogen infection of plants [63]. In cultured potato cells treated with an elicitor that causes late blight disease, there was an accumulation of 9-LOX transcripts, which indicated plant defences [64]. In our analysis, LOC107829039 and LOC107770253 are members of the lipoxygenase family. In addition, LOC107774644 might encode an isoform of the MADS-box TFs involved in stress responses, and a MADS-box TF in Capsicum annuum was induced by MeJA [64]. LOC107829039, LOC107770253, and LOC107774644 all have cis-acting elements upstream of their respective lncRNAs associated with plant defence and stress ( Figure 6). Therefore, we hypothesized that MeJA was produced after plant stress and that lncRNAs may be activated by their upstream MeJA-responsive cis-acting elements. Subsequently, the presence of plant defence-related elements allows lncRNAs and their target genes to engage in plant disease resistance. Moreover, many genes encoding the transcription factors have also become targets of MeJA-responsive lncRNAs, such as RAX3-like (LOC107785742), RAX1-like (LOC107788895), myb-related protein 306-like (LOC107775566), and MADS-box 23-like isoform X1 (LOC107774644). In summary, these TFs are highly likely to regulate MeJA-responsive lncRNAs and their target genes at the transcriptional level.

Plant Materials and Growth Conditions
Suspension cultures of Nicotiana tabacum L. cv. Bright Yellow-2 (BY-2) cells [65] were grown and maintained, according to that which was previously reported [57]. Briefly, BY-2 cells were cultured in a 50-mL half-strength Murashige and Skoog (MS) medium (pH 5.8), supplied with 3% (w/v) sucrose and 0.2 mg/L of 2,4-dichlorophenyoxyacetic acid (2,4-D) in a 250-mL flask, shook with 150 rpm at 28 • C, and subcultured in a fresh medium for every week. For the MeJA treatment, a four-day cell suspension was subcultured in a fresh medium without 2,4-D for one day, then treated with MeJA at a final concentration of 100 µM [dissolved in dimethyl sulfoxide (DMSO) and sterilised with a 0.22-µm Millex syringe filter] for 2 h. Cells treated with an equal volume of DMSO were used as a mock treatment. Vacuum filtration was used to collect cells from the treated cultures for further analysis.

General RNA-Seq Read Mapping and Assembly
The paired-end transcriptome sequencing data of BY-2 cells were obtained from the NCBI Sequence Read Archive (SRA: SRP026451). The general RNA-seq processing pipeline is shown in the workflow ( Figure 1A). We used FastQC (http://www.bioinformatics. babraham.ac.uk/projects/download.html#fastqc) (version 0.11.9) (the Babraham Institute, Cambridge, UK) and a FastX-Toolkit (hannonlab.cshl.edu/fastx toolkit/) (version 0.0.14) (Cold Spring Harbor Laboratory, Laurel Hollow, NY, USA) to assess the quality of the initial reads and trimmed the raw data. The minimal quality score to preserve was set to 20, the minimum percent of the bases was set to 90, and the first seven nucleotides were trimmed off. For indexing, the reference genome sequences of N. tabacum (GCF 000715135.1) from Bowtie2 (http://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.6/) (version 2.2.6) (Johns Hopkins University, Baltimore, MD, USA) were used [66]. TopHat2 (https: //ccb.jhu.edu/software/tophat/index.shtml) (version 2.0.14) (University of Maryland, College Park, MD, USA) aligned clean reads to the genome and discovered transcript splice sites (maximum intron size = 5000) [67]. Cufflinks (http://cole-trapnell-lab.github. io/cufflinks/install/) (version 2.2.1) (University of California, Berkeley, CA, USA) were utilized for the transcript assembly, combination, and comparison, in which the transcript expression levels were normalized using fragments per kilobase of transcript per million fragments (FPKM).

Prediction of Cis-Regulating Target Genes of MeJA-Responsive lncRNAs
Based on the cis-acting mechanism, we predicted the target genes of all lncRNAs. Potential cis-regulating target genes were those transcribed within a 100 kb window upstream or downstream of lncRNAs [68]. Bedtools (https://github.com/arq5x/bedtools2/releases) (version 2.30.0) (The University of Utah, Salt Lake City, UT, USA) was used to find the lncRNAs' nearby protein-coding genes (PCGs). Venns were utilized to show up-and downregulated MeJA-responsive lncRNAs, cis-regulating lncRNAs, and their target genes.

Expression Pattern Comparison of Cis-Regulated Target Genes with Their MeJA-Responsive lncRNAs
We analysed all upregulated DEgenes, all downregulated DEgenes, all cis-regulated target genes of upregulated MeJA-responsive lncRNAs, and all cis-regulated target genes of downregulated MeJA-responsive lncRNAs using Venn modules from the TBtools software.

Gene Ontology and Pathway Enrichment Analysis
We concentrated on the potential biological functions of cis-regulated target genes of MeJA-responsive lncRNAs (including both down-and upregulated target genes). Using an eggNOG-mapper (http://eggnog-mapper.embl.de/) (EMBL), we performed a Gene Ontology (GO) analysis of the protein sequences encoded by these target genes. Using the KEGG (Kyoto Encyclopedia of Genes and Genomes) database (http://www.genome.ad.jp/ kegg/) (Kyoto University, Kyoto, Japan), enrichment pathway analysis was performed to investigate the probable activities of the target genes in plant pathways.

Cis-Regulatory Element Prediction in the Promoter Sequences
To understand the cis-regulatory element status of the MeJA-responsive lncRNAs, we implemented cis-regulatory element analysis with PlantCARE (http://bioinformatics.psb. ugent.be/webtools/plantcare/html/) (Universiteit Gent, Ghent, Belgium). The 2000 bp upstream of the lncRNA gene sequences represented the promoter region.

Conclusions
In this study, we identified 7848 Nicotiana genus-specific lncRNAs in BY-2 cells, of which 629 differentially expressed (DE) lncRNAs were defined as MeJA-responsive lncRNAs. The combined analysis of the cis-regulated lncRNAs and their target genes revealed potential up-and downregulated target functional genes. Additionally, a series of MeJA-responsive lncRNA target genes played a role in nicotine synthesis and disease and insect resistance, which was associated with their MeJA-and defence-related cis-regulatory elements. In conclusion, the foregoing outcome indicates that MeJA-responsive lncRNAs regulate nicotine biosynthesis and disease resistance by regulating their potential target genes in BY-2 cells. While these findings reveal a wide range of MeJA-responsive lncRNAs in tobacco BY-2 cells, further research is required to elucidate the molecular mechanisms of coding and noncoding genes at the transcription level and protein level and in epigenetic inheritance. In brief, this work provided novel insights into their plant defence-associated target genes and an abundant gene resource in tobacco for functional research.

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