Systematic Identiﬁcation of Long Non-Coding RNAs under Allelopathic Interference of Para-Hydroxybenzoic Acid in S. lycopersicum

: The importance of long noncoding RNAs (lncRNAs) in plant development has been established, but a systematic analysis of the lncRNAs expressed during plant allelopathy has not been carried out. We performed RNA-seq experiments on S. lycopersicum subjected to different levels of para-hydroxybenzoic acid (PHBA) stress during plant allelopathy and identiﬁed 61,729 putative lncRNAs. Of these, 7765 lncRNAs cis-regulated 5314 protein-coding genes (PGs). Among these genes, 1116 lncRNAs and 2239 PGs were involved in a complex web of transcriptome regulation, and we divided these genes into 12 modules. Within these modules, 458 lncRNAs and 975 target genes were found to be highly correlated. Additionally, 989 lncRNAs trans-regulated 1765 PGs, and we classiﬁed them into 11 modules, within which 335 lncRNAs were highly correlated with their 633 corresponding target genes. Only 98 lncRNAs in S. lycopersicum had homologs in the lncRNA database of Arabidopsis thaliana , all of which were affected by the PHBA treatments. MiRNAs that interacted with both mRNAs and lncRNAs were selected on the basis of weighted correlation network analysis (WGCNA) results to make lncRNA-miRNA-mRNA triplets. Our study presents a systematic identiﬁcation of lncRNAs involved in plant allelopathy in S. lycopersicum and provides research references for future studies.


Introduction
Since the proposal of the central dogma of molecular biology in 1961, RNA has been considered an intermediate that translates genetic information from DNA to protein [1,2]. However, with the discovery of noncoding RNAs (ncRNAs), we found that intermediate RNAs represent a fraction of all RNAs [2,3]. As the main products of the eukaryotic transcriptome, ncRNAs function as structural, catalytic, or regulatory RNAs rather than as protein encoders [2,[4][5][6]. Based on RNA length, regulatory ncRNAs can be further classified as small ncRNAs (<200 bp, e.g., miRNAs, siRNAs, and piRNAs) or long ncRNAs (lncRNAs) (>200 bp, e.g., lincRNAs and macroRNAs) [2,7]. In the early 1990s, with the appearance of the X-chromosome-silencing phenomenon, long ncRNAs (lncRNAs) were first discovered, having once been considered the 'dark matter' of transcriptomes [8]. Recently, lncRNAs have been discovered to have strong and universal regulatory effects on gene expression at the post-transcriptional, transcriptional, and epigenetic levels [9][10][11][12][13].
lncRNAs play important roles in many biological processes in plants, especially in developmental regulation and stress responses [14][15][16]. A total of 1832 lncRNAs in Arabidopsis were changed after cold, drought, high-salt, and abscisic acid (ABA) treatments [17]. Several lncRNAs, such as COLDAIR and COOLAIR, positively respond to cold and function in flowering induction [18,19]. IPS1 and At4, induced by phosphate starvation, regulate the shoot dynamic balance of phosphate through the blockage of the repression of miR399 on its target gene, PHO2 [20][21][22][23]. Further exploration of the function of lncRNAs in plant responses to stress would provide approaches to discover the stress response networks [24].
As an important factor affecting agricultural production, allelopathy has become a hot topic in ecology, horticulture, agronomy, and other research areas [25][26][27][28][29]. Allelochemicals have been proven to be responsible for numerous biochemical and physiological changes that cause plant allelopathy [30,31]. The plant growth regulatory system, respiration, photosynthesis, the antioxidant system, and water and nutrient uptake are the key physiological and biochemical systems and processes in plants in which changes are induced by allelochemicals [32,33].
Transcription factors, DNA methylation events, chromatin modifications, and mi-croRNAs have been discovered in the gene regulatory processes described above [34][35][36][37]. LncRNAs are critical in plant gene expression regulation in response to stress [38][39][40], which indicates that lncRNAs might also function in plant response to allelopathy. The involvement of lncRNAs in gene transcription regulation progress have not been explored. The expression of most lncRNAs is tissue-specific, which allows them to be discovered by transcriptome sequencing (RNA-seq). By examining RNA-seq data of leaves under salt stress, long noncoding RNAs were identified and characterized in Medicago truncatula [41].
However, plant allelopathy is a complex process. To explore the allelopathy network in plants, it is necessary to study the dynamic changes in gene expression under different degrees of stress. The roles of lncRNAs in plant allelopathy in S. lycopersicum were explored in this study. Firstly, lncRNAs expressed in S. lycopersicum leaves under different PHBA treatments were identified with RNA-seq. Next, cis-and trans-regulated target PGs (TTPGs) for lncRNAs were predicted, and the coexpressions between lncRNAs and their targets were analyzed. Third, the conservation of lncRNAs in S. lycopersicum and Ahrabidopsis thaliana was compared. Lastly, we explored the possibility that these lncRNAs might be endogenous pseudotarget mimics (eTMs) of known S. lycopersicum miRNAs. Overall, this study helps to understand the role of lncRNAs in the molecular network of plant allelopathy and gives new clues about plant response to stress.

Plant Materials and Growth Conditions
The S. lycopersicum cultivar 'Diana' was grown at the experimental farm of the Shandong Facility Horticultural Laboratory, Weifang, Shandong Province, China. Tomato seeds were planted in plug plates (Wei Nong Company, Taizhou, China). One month later, we treated the S. lycopersicum plants with different concentrations of PHBA (0 mmol/L, 15 mmol/L, and 30 mmol/L) for 2 days. Then, the leaves from the middles of the plants were collected as samples. Specifically, the samples for the CK, T1, and T2 treatments corresponded to the 0 mmol/L, 15 mmol/L, and 30 mmol/L concentrations of PHBA, respectively ( Figure 1A). We collected leaf samples from 60 individual tomato seedlings and stored them at −75 • C in an ultra-low temperature freezer for the following experiments.

RNA Sequencing and Identification of lncRNAs
Total RNA was extracted from the leaf samples for the different PHBA treatments described above using the MagMAX Plant RNA Isolation Kit (Thermo Fisher Scientific, Carlsbad, CA, USA). RNA quality was measured with a Thermo Scientific NanoDrop 2000 system (Thermo Fisher Scientific, USA), and RNA integrity was examined with an Agilent 2100 Bioanalyzer System (Agilent Technologies, Santa Clara, CA, USA) [42]. RNAs with integrity numbers greater than 7 were selected for strand-specific RNA-seq library construction. First, ribosomal RNAs were removed from the total RNAs using the plant Ribo-Zero rRNA Removal Kit (Illumina, San Diego, CA, USA). A strand-specific RNA-seq library was constructed with the Ribo-Zero RNA fraction via the Illumina TruSeq RNA Library Preparation Kit v2 and sequencing on an Illumina HiSeq 4000 instrument.

RNA Sequencing and Identification of lncRNAs
Total RNA was extracted from the leaf samples for the different PHBA treatments described above using the MagMAX Plant RNA Isolation Kit (Thermo Fisher Scientific, Carlsbad, CA, USA). RNA quality was measured with a Thermo Scientific NanoDrop 2000 system (Thermo Fisher Scientific, USA), and RNA integrity was examined with an Agilent 2100 Bioanalyzer System (Agilent Technologies, Santa Clara, CA, USA) [42]. RNAs with integrity numbers greater than 7 were selected for strand-specific RNA-seq library construction. First, ribosomal RNAs were removed from the total RNAs using the plant Ribo-Zero rRNA Removal Kit (Illumina, San Diego, CA, USA). A strand-specific RNA-seq library was constructed with the Ribo-Zero RNA fraction via the Illumina TruSeq RNA Library Preparation Kit v2 and sequencing on an Illumina HiSeq 4000 instrument.
After removing adaptors and low-quality reads from raw reads using Trimmomatic [43], 1,003,711,588 paired-end clean reads were produced (111,523,510 on average for each sample). They were subjected to de novo transcript assembly using Trinity (v2.8.4, Manfred G. Grabherr, Cambridge MA, USA), with strand-specific RNA-Seq read orientation [44]. Duplicates were then removed from the assembled transcripts using CD-HIT (v4.7, Limin Fu, La Jolla CA, USA) with default parameters [45]. Finally, we obtained 82,126 unigenes longer than 200 bp for lncRNA screening. The lncRNA identification process is shown in Figure S1. To improve the accuracy of lncRNA detection, three software programs, RNAplonc (v1.0, Tatianne da Costa Negri, São Paulo, Brazil) [46], PLncPRO After removing adaptors and low-quality reads from raw reads using Trimmomatic [43], 1,003,711,588 paired-end clean reads were produced (111,523,510 on average for each sample). They were subjected to de novo transcript assembly using Trinity (v2.8.4, Manfred G. Grabherr, Cambridge MA, USA), with strand-specific RNA-Seq read orientation [44]. Duplicates were then removed from the assembled transcripts using CD-HIT (v4.7, Limin Fu, La Jolla CA, USA) with default parameters [45]. Finally, we obtained 82,126 unigenes longer than 200 bp for lncRNA screening. The lncRNA identification process is shown in Figure S1. To improve the accuracy of lncRNA detection, three software programs, RNAplonc (v1.0, Tatianne da Costa Negri, São Paulo, Brazil) [46], PLncPRO (v1.1, Urminder Singh, New Delhi, India) [47], and FEELnc (v0.1.1, Valentin Wucher, Rennes, Cedex, France) [48], with default parameters, were employed to identify lncRNAs. The lncRNAs obtained with these three tools were compared to obtain the shared pre-lncRNA candidates. To further rule out the possibility that the pre-lncRNAs were from coding genes, these pre-lncRNAs were blasted against coding DNA sequences of S. lycopersicum (SL.3.0) using BLAT v36 × 2 software [49] with an E-value of 10 −5 . The filtered lncRNA candidates were then localized to the S. lycopersicum reference genome using BLAT to remove repeat sequences, so we could eliminate ones that had more than nine positions and which were considered to be repeats in the genome. Eventually, we identified highly reliable putative lncRNAs in S. lycopersicum.

Expression Analysis of Putative lncRNAs and Quantitative Reverse Transcription PCR (qRT-PCR) Validation
RSEM (v1.3.3, Bo Li, Madison, WI, USA) was used to compare the sequenced reads with the transcripts to calculate the gene expression amounts and carry out a differential gene expression analysis. The levels of gene expression were reflected by FPKM values, and a threshold of FPKM > 0.1 was used for the expressed genes. The DESeq package in R software (v3.6.1, R Core Team, Vienna, Austria) was used for pairwise differential expression analysis. The website http://bioinformatics.psb.ugent.be/webtools/Venn/ (accessed on 14 June 2019) was used for Venn diagram drawing [50]. qRT-PCR was performed to validate the expressions of the putative lncRNAs with the same RNA samples as those used for the RNA-seq [51]. The primers are listed in Table S1. GAPDH was chosen as the reference gene. The DDCt method was used to calculate the results [52].

Target Gene Prediction and Coexpression Analysis of lncRNAs with PGs
WGCNA was used to explore the cis-regulatory relationships between lncRNAs and PGs by calculating the correlations with chromosomal neighboring genes within 100 kb from the middle of a given lncRNA [50,53]. RNAplex software was used to predict the trans-acting lncRNAs and their target PGs (TPGs) [54]. Their coexpression was analyzed using WGCNA (cis-analysis: soft threshold = 21; trans-analysis: soft threshold = 22). R was used to calculate the Spearman correlations for the expression of lncRNAs and their target gene pairs.
The GO Analysis Toolkit and the Database for the Agricultural Community (agriGO; http://bioinfo.cau.edu.cn/agriGO/analysis.php) (accessed on 19 July 2019) were used for the GO analysis of the TPGs. The hypergeometric test was used as the statistical method, while Yekutieli (FDR under dependency) was employed as the multi-test adjustment method [50].

Similarity Analysis of lncRNAs from S. lycopersicum and Arabidopsis thaliana
The Plant Long noncoding RNA Database (PLncDB, http://chualab.rockefeller.edu/ gbrowse2/homepage.html) (accessed on 16 July 2019) was used to download the lncRNA data for Arabidopsis thaliana [50]. BLAST software (E-value set to 10 −5 ) was used to compare the conservation and similarity levels of lncRNAs between S. lycopersicum and Arabidopsis thaliana [50].

Prediction of Endogenous Target Mimics for miRNAs
MicroRNA data for S. lycopersicum were obtained from miRbase (http://www.mirbase. org/) (accessed on 26 July 2019). WGCNA was used to predict the relationships between mRNAs, lncRNAs, and miRNAs. A threshold weight value of >0.15 was used for the selection of lncRNA-miRNA-mRNA triplets. Cytoscape (https://cytoscape.org/) (accessed on 26 July 2019) was used to plot the networks of lncRNAs, miRNAs, and mRNAs.

Identification of lncRNAs in S. lycopersicum
Substantial numbers of novel transcripts could be discovered via the combination of high-throughput RNA-seq and de novo assembly. In this study, the 'Diana' cultivar of S. lycopersicum was used for RNA-seq. As a new cultivar of S. lycopersicum, 'Diana' is widely cultivated in China. S. lycopersicum; an important crop species, it is also affected by continuous cropping obstacles, which are mainly caused by allelopathy. PHBA, a common allelochemical, was studied in depth, and the results showed that photosynthesis, respiratory effects, and gene expression in the leaves were affected by PHBA. Thus, the 'Diana' cultivar of S. lycopersicum was chosen to study gene expression in PHBA-induced allelopathy. We built the unigene library with three samples, including leaves undergoing three different treatments (CK: leaves of S. lycopersicum under 0 mmol/L PHBA treatment for 2 days; T1: leaves of S. lycopersicum under 15 mmol/L PHBA treatment for 2 days; and T2: leaves of S. lycopersicum under 30 mmol/L PHBA treatment for 2 days) ( Figure 1A).
A total of 82,126 unigenes with an overall length greater than 200 bp and poly(A) tails were identified after de novo assembly. Following the protocol represented in Figure S1, 65,178 unigenes were acquired as credible lncRNAs, among which 61,729 lncRNAs were expressed during the plant allelopathy progression in S. lycopersicum in response to PHBA (Table S2). Ten lncRNAs were randomly chosen, and their expression was verified by qRT-PCR; all results were consensual ( Figure S2).
The average length of the 61,729 assumed lncRNAs was 565 bp, the minimum length was 201 bp, and the maximum length was 19,161 bp ( Figure 1B and Table S2). Of these lncRNAs, 80.68% were shorter than 600 bp, while 9.59% were longer than 1000 bp. About two-thirds of them expressed at low levels, with FPKM values (log2) ranging from 0 to 3 ( Figure 1C and Table S2). Interestingly, the trends in the data for the lncRNAs corresponded to those for the PGs across all of the samples ( Figure 1D and Table S3).

Expression of lncRNAs in S. lycopersicum Leaves under Different PHBA Treatments
As we knew, 61,729 lncRNAs were expressed in S. lycopersicum leaves during the different PHBA treatments, among which 34,000 lncRNAs were expressed throughout all three treatments. Almost half of the lncRNAs were expressed exclusively in two stages, while a large number of lncRNAs were expressed in only one stage (Figure 2A). We compared the differences in expression between the lncRNAs and PGs in S. lycopersicum leaves under different PHBA treatments (Table S4). The numbers of differentially expressed PGs and lncRNAs between the treated and untreated S. lycopersicum leaves increased gradually with increasing PHBA concentration and reached an apex under the T2 treatment ( Figure 2B). We compared the differences in expression between the lncRNAs and PGs in S. lycopersicum leaves under different PHBA treatments (Table S4). The numbers of differentially expressed PGs and lncRNAs between the treated and untreated S. lycopersicum leaves increased gradually with increasing PHBA concentration and reached an apex under the T2 treatment ( Figure 2B).

Prediction of Cis-Regulated Target PGs (CTPGs) of Lncrnas
The expression of proximal and distal PGs has been found to be regulated by lncRNAs through cis-and trans-acting mechanisms [55]. The regulation of genes located on the same chromosome as lncRNAs can be ascribed to cis-regulation [56,57], the investigation of which relies on well-annotated genomes.
The proximal PGs located in a 100 kb genomic window were searched as cis-regulated target genes of lncRNAs [53]. A total of 7765 lncRNAs were found to have potential cisregulatory effects on 5314 PGs in 38,845 gene pairs altogether (Table S5), among which, 89% of lncRNAs targeted more than one gene, 71% targeted two to seven genes, while only two lncRNAs targeted fourteen genes ( Figure 3A). Approximately 5% of the PGs corresponded to only one lncRNA, and only six PGs were cis-regulated by up to nineteen lncRNAs ( Figure 3A).

Prediction of Cis-Regulated Target PGs (CTPGs) of Lncrnas
The expression of proximal and distal PGs has been found to be regulated by lncRNAs through cis-and trans-acting mechanisms [55]. The regulation of genes located on the same chromosome as lncRNAs can be ascribed to cis-regulation [56,57], the investigation of which relies on well-annotated genomes.
The proximal PGs located in a 100 kb genomic window were searched as cis-regulated target genes of lncRNAs [53]. A total of 7765 lncRNAs were found to have potential cis-regulatory effects on 5314 PGs in 38,845 gene pairs altogether (Table S5), among which, 89% of lncRNAs targeted more than one gene, 71% targeted two to seven genes, while only two lncRNAs targeted fourteen genes ( Figure 3A). Approximately 5% of the PGs corresponded to only one lncRNA, and only six PGs were cis-regulated by up to nineteen lncRNAs ( Figure 3A). Subsequently, the coexpression relationships between lncRNAs and protein-coding gene pairs were studied using WGCNA. A total of 1116 lncRNAs and 2239 PGs were involved in transcriptome regulation, which formed a complex network (Data S1), including 12 modules ( Figure 4A and Table S6). Different modules contained different proportions of lncRNAs, ranging from 16.7% in Module 10 (M10) to 57.6% in M6, with an average of 38.2%. Gene expression patterns for each module were distinguished ( Figure 5). Signifi- Subsequently, the coexpression relationships between lncRNAs and protein-coding gene pairs were studied using WGCNA. A total of 1116 lncRNAs and 2239 PGs were involved in transcriptome regulation, which formed a complex network (Data S1), including 12 modules ( Figure 4A and Table S6). Different modules contained different proportions of lncRNAs, ranging from 16.7% in Module 10 (M10) to 57.6% in M6, with an average of 38.2%. Gene expression patterns for each module were distinguished ( Figure 5). Significant GO terms were identified in each module, ranging from 11 in M11 to 133 in M1 (Table S7). GO:0032957, which is associated with the inositol trisphosphate metabolic process, was found in M7 and M11, while GO:0050789 (regulation of biological process) was found in M1 and M6. Among these modules, the highly expressed genes in M2 were involved in photosynthesis and cell wall organization and were annotated as being associated with allelopathy ( Figure 6A). The genes in M13 were also predicted to be involved in photosynthesis and cell wall biogenesis ( Figure 6B). S7). GO:0032957, which is associated with the inositol trisphosphate metabolic process, was found in M7 and M11, while GO:0050789 (regulation of biological process) was found in M1 and M6. Among these modules, the highly expressed genes in M2 were involved in photosynthesis and cell wall organization and were annotated as being associated with allelopathy ( Figure 6A). The genes in M13 were also predicted to be involved in photosynthesis and cell wall biogenesis ( Figure 6B).   S7). GO:0032957, which is associated with the inositol trisphosphate metabolic process, was found in M7 and M11, while GO:0050789 (regulation of biological process) was found in M1 and M6. Among these modules, the highly expressed genes in M2 were involved in photosynthesis and cell wall organization and were annotated as being associated with allelopathy ( Figure 6A). The genes in M13 were also predicted to be involved in photosynthesis and cell wall biogenesis ( Figure 6B).     In addition, we calculated the Spearman correlations (correlation coefficient cutoff = 0.9) for the expression of each lncRNA and its CTPG pair. The expressions of 458 lncRNAs and their corresponding 975 CTPGs were found to be strongly correlated (Table S8). A total of 1092 gene pairs were positively correlated, and 963 pairs were negatively correlated. Among these target genes, 1348 were upregulated in S. lycopersicum leaves in at least one treatment, and 892 were downregulated (Table S9).
Using WGCNA, 989 lncRNAs and their 1765 TTPGs were divided into 11 modules ( Figure 4B). The proportions of lncRNAs ranged from 11.1% in M14 to 67.1% in M15, with an average of 35.9% in these modules (Table S6). Except for M15, M19, M20, M21, and M23, remarkable GO terms were verified in most modules, with a range from 106 in M15 to 6 in M23 (Table S11). Distinguishable expression patterns were found in all of the modules (Figure 7). For example, genes in M16 and M19 were extremely highly expressed in plants undergoing the T2 treatment ( Figure 7) and were annotated as being associated with the following terms: reactive oxygen species metabolic process, antibiotic catabolic process, and hydrogen peroxide catabolic process (Figure 8). Regulatory networks that were widely interconnected between mRNAs and lncRNAs were found in all modules of trans-regulatory activity (Data S1).
In addition, we calculated the Spearman correlations (correlation coefficient cutoff = 0.9) for the expression of each lncRNA and its CTPG pair. The expressions of 458 lncRNAs and their corresponding 975 CTPGs were found to be strongly correlated (Table S8). A total of 1092 gene pairs were positively correlated, and 963 pairs were negatively correlated. Among these target genes, 1348 were upregulated in S. lycopersicum leaves in at least one treatment, and 892 were downregulated (Table S9).
Using WGCNA, 989 lncRNAs and their 1765 TTPGs were divided into 11 modules ( Figure 4B). The proportions of lncRNAs ranged from 11.1% in M14 to 67.1% in M15, with an average of 35.9% in these modules (Table S6). Except for M15, M19, M20, M21, and M23, remarkable GO terms were verified in most modules, with a range from 106 in M15 to 6 in M23 (Table S11). Distinguishable expression patterns were found in all of the modules (Figure 7). For example, genes in M16 and M19 were extremely highly expressed in plants undergoing the T2 treatment ( Figure 7) and were annotated as being associated with the following terms: reactive oxygen species metabolic process, antibiotic catabolic process, and hydrogen peroxide catabolic process (Figure 8). Regulatory networks that were widely interconnected between mRNAs and lncRNAs were found in all modules of trans-regulatory activity (Data S1).   A total of 335 trans-acting lncRNAs were highly correlated with their 633 ta (Spearman coefficient rho (rs) > 0.9) (Table S12); 681 gene pairs were positively correl while 554 pairs were negatively. The expression patterns of the targets indicated that were upregulated and that 730 were downregulated in at least one treatment (Table A total of 335 trans-acting lncRNAs were highly correlated with their 633 targets (Spearman coefficient rho (rs) > 0.9) (Table S12); 681 gene pairs were positively correlated, while 554 pairs were negatively. The expression patterns of the targets indicated that 1035 were upregulated and that 730 were downregulated in at least one treatment (Table S13).

Similarity Alignment and Conservation Analysis of lncRNAs in S. lycopersicum and Arabidopsis thaliana
As a model plant and close species of S. lycopersicum, Arabidopsis thaliana was used to analyze the lncRNA similarities between them. Although there were 65,178 lncRNAs in S. lycopersicum, only 98 had homologs in Arabidopsis thaliana (Table S14). Expression analysis revealed that these 98 lncRNAs were affected by the PHBA treatments (Figure 9).

Similarity Alignment and Conservation Analysis of lncRNAs in S. lycopersicum and Arabidopsis thaliana
As a model plant and close species of S. lycopersicum, Arabidopsis thaliana was used to analyze the lncRNA similarities between them. Although there were 65,178 lncRNAs in S. lycopersicum, only 98 had homologs in Arabidopsis thaliana (Table S14). Expression analysis revealed that these 98 lncRNAs were affected by the PHBA treatments ( Figure 9).

Prediction of lncRNAs as Endogenous Target Mimics of miRNAs
Recent research has indicated that lncRNAs interact with miRNAs as competitive endogenous RNAs (ceRNAs) and participate in expression regulation of targets [58,59]. The potential lncRNAs were predicted to be eTM miRNAs using S. lycopersicum miRNAs which were previously obtained from miRbase (http://www.mirbase.org/, accessed on 26 July 2019).
The miRNAs that interacted with both mRNAs and lncRNAs were selected from the WGCNA results to make lncRNA-miRNA-mRNA triplets (Table S15). A weight value greater than 0.15 was used to select the lncRNA-miRNA-mRNA triplets (Table S16). Cytoscape was used to plot the interactions between lncRNAs, miRNAs, and mRNAs ( Figure  10).

Prediction of lncRNAs as Endogenous Target Mimics of miRNAs
Recent research has indicated that lncRNAs interact with miRNAs as competitive endogenous RNAs (ceRNAs) and participate in expression regulation of targets [58,59]. The potential lncRNAs were predicted to be eTM miRNAs using S. lycopersicum miRNAs which were previously obtained from miRbase (http://www.mirbase.org/, accessed on 26 July 2019).
The miRNAs that interacted with both mRNAs and lncRNAs were selected from the WGCNA results to make lncRNA-miRNA-mRNA triplets (Table S15). A weight value greater than 0.15 was used to select the lncRNA-miRNA-mRNA triplets (Table S16). Cytoscape was used to plot the interactions between lncRNAs, miRNAs, and mRNAs ( Figure 10).

Discussion
Although genome-wide lncRNA searches and related studies have been carried out in Arabidopsis thaliana [60], rice [61], maize [62], wheat [63], and Populus trichocarpa , identification of lncRNA functions in plant response to allelopathy is lacking. In this study, an RNA-seq experiment was carried out on tomato plants treated with different concentrations of PHBA, which systematically identified plant allelopathy-related lncRNAs. The tomato material we used in this study is a new variety that is widely planted in North China and thus has an important reference value for plant allelopathy in vegetable production. In addition, the differentially expressed lncRNAs under different PHBA treatments were valuable for studying plant allelopathy in S. lycopersicum. Moreover, as an important vegetable crop grown worldwide, S. lycopersicum is considered a reference species for Solanum genomic studies [64]. The genome-wide characterization of lncRNAs in S. lycopersicum will provide a unique annotation resource for the S. lycopersicum genome.
Gene expression can be cis-or trans-regulated by lncRNAs [65]. Cis-acting lncRNAs control gene expression around their transcription sites, while trans-acting lncRNAs regulate the expression of genes at independent loci [66][67][68]. Additionally, the overlap of some cis-and trans-acting lncRNAs increases the complexity of the molecular roles of

Discussion
Although genome-wide lncRNA searches and related studies have been carried out in Arabidopsis thaliana [60], rice [61], maize [62], wheat [63], and Populus trichocarpa, identification of lncRNA functions in plant response to allelopathy is lacking. In this study, an RNA-seq experiment was carried out on tomato plants treated with different concentrations of PHBA, which systematically identified plant allelopathy-related lncRNAs. The tomato material we used in this study is a new variety that is widely planted in North China and thus has an important reference value for plant allelopathy in vegetable production. In addition, the differentially expressed lncRNAs under different PHBA treatments were valuable for studying plant allelopathy in S. lycopersicum. Moreover, as an important vegetable crop grown worldwide, S. lycopersicum is considered a reference species for Solanum genomic studies [64]. The genome-wide characterization of lncRNAs in S. lycopersicum will provide a unique annotation resource for the S. lycopersicum genome.
Gene expression can be cis-or trans-regulated by lncRNAs [65]. Cis-acting lncRNAs control gene expression around their transcription sites, while trans-acting lncRNAs regulate the expression of genes at independent loci [66][67][68]. Additionally, the overlap of some cis-and trans-acting lncRNAs increases the complexity of the molecular roles of lncR-NAs [50,69]. Since lncRNAs are often coexpressed with their TPGs, a method for inferring the putative biological functions of lncRNAs was developed based on the relationships between the expressions of lncRNAs and TPGs [50,70]. The CTPGs and TTPGs of lncRNAs were predicted, and a coexpression analysis was performed [50]. Strong coexpression between cis-or trans-acting lncRNAs and their TPGs was found in this study. To date, no lncRNA has been characterized in plant allelopathy processes. Further investigation of these lncRNAs that are highly coexpressed with TPGs with known functions may provide some clues that will enable discovery of the functions of lncRNAs in gene expression regulation of allelopathy.
Compared with the study of functionality, the study of evolution needs to be further explored. For plant lncRNAs, a recent BLAST search of B. rapa lncRNAs against the lncRNAs of Arabidopsis thaliana yielded few homologs (approximately 5%) [50], indicating that evolutionary conservation is limited in plant lncRNAs [50]. In this study, we found that only approximately 0.15% (98 out of 65,178) of lncRNAs in S. lycopersicum had homologs in Arabidopsis thaliana. This result also suggested limited evolutionary conservation in plant lncRNAs.
Recent studies have found lncRNAs involved in miRNA regulation mechanisms as ceRNAs, which were initially studied in Arabidopsis thaliana and tobacco [71,72]. In this study, we predicted the relationships between lncRNAs, miRNAs, and mRNAs using the miRNAs that were previously identified in S. lycopersicum. We propose that certain interactions between these lncRNAs, mRNAs, and miRNAs play fundamental roles in plant allelopathy.

Conclusions
In conclusion, screening and functional analysis of the genome of S. lycopersicum enabled the identification of a set of lncRNAs. The analysis of their functional links with mRNAs and miRNAs and the similarity alignment and conservation analyses of lncRNAs in S. lycopersicum and Arabidopsis thaliana may provide some molecular clues about plant allelopathy due to PHBA in S. lycopersicum.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/horticulturae8121134/s1, Data S1: The coexpression relationships between lncRNAs and protein-coding gene pairs, Figure S1: Informatics pipeline for the identification of lncRNAs in S. lycopersicum, Figure S2: Experimental validation of the lncRNA expression by qRT-PCR. (The horizontal axis represents different treatments of PHBA, the vertical axis represents the relative expression level of lncRNAs.), Table S1: The primers used in this study, Table S2: LncRNAs identified in S. lycopersicum, Table S3: Numbers of lncRNAs and PGs detected in S. lycopersicum, Table S4: Differential expressions of lncRNAs and mRNAs, Table S5: LncRNAs and their CTPGs,  Table S6: The weighted gene co-expression network analysis modules of lncRNAs and corresponding TPGs, Table S7: GO terms and annotations of cis-regulated modules, Table S8: LncRNAs and cis-regulated target mRNAs that have correlation coefficients greater than 0.9 and p-values less than 1.0 × 10 −5 , Table S9: Cis-regulated target mRNAs that are differentially expressed, Table S10: LncRNAs and their trans-regulated corresponding target mRNAs analyzed using LncTar, Table S11: GO terms and annotations of trans-regulated modules, Table S12: LncRNAs and trans-regulated target mRNAs that have correlation coefficients greater than 0.9 and p-values less than 1.0 × 10 −5 , Table S13: Trans-regulated target mRNAs that are differentially expressed, Table S14: Similarity alignment of lncRNAs in S. lycopersicum and Arabidopsis thaliana, Table S15: MiRNA target prediction by psRNATarget, Table S16: The lncRNA-miRNA-mRNA pairs with weight values greater than 0.15.