The Gossypium hirsutum genome sequences and the proteome sequences were downloaded from the CottonGen database (http://www.cottongen.org) (Yu et al., 2014). The amino acid sequences of TOR from A. thaliana, which were acquired from The Arabidopsis Information Resource version 10 (TAIR 10) (http://www.arabidopsis.org), were used as query sequences to identify complete TOR members in the G. hirsutum protein database through the blastp program. The full-length protein sequences of Arabidopsis and G. hirsutum TOR were aligned with Clustal W, and MEGA 5.0 software was applied to construct a neighbor-joining phylogenetic tree using the bootstrap method with 1,000 replicates (Tamura et al., 2011). The positions of the exons and introns were acquired from the gff3 file using a Perl script, and the online tool GSDS 2.0 was used to display the gene structure (Hu et al., 2015).

In this study, upland cotton (Gossypium hirsutum L. cv CCRI24) was used. The cotton seeds were surface-sterilized according to the method described previously (Yang et al., 2014). The sterilized seeds were germinated on a Murashige and Skoog (MS) medium (with 1% agar) supplemented with DMSO (as a control) or 5 μM AZD (Selleckchem, Houston, TX, USA). Then they were grown in a controlled environment at 25°C under 16 h continuous light and 8 h darkness for 4, 6, and 8 days, respectively. These cotton seedlings were photographed. The primary root lengths, fresh weights and lateral root numbers were measured manually. Three biological experiments were performed, and each consisted of 20 plants per treatment. The data was expressed as the mean ± SD of three independent experiments.

The conserved cDNA fragments of GhTOR genes were amplified and cloned into the pTRV2 vectors according to the method described in the previous study (Liu et al., 2002). The recombinant plasmids were transformed into A. tumefaciens strain GV3101. Cotton VIGS experiments were performed following the previously described procedure (Pang et al., 2013). Cotyledons of 7-day-old CCRI24 cotton seedlings were used for the injections. All plants infiltrated with the empty pTRV2 were used as the negative controls in this research. Plants infiltrated with the pTRV2 vector containing the cDNA fragment of GhPDS gene were used as the positive controls to test the efficiency of VIGS, and the albino phenotype was examined (Pang et al., 2013). Three biological experiments were performed, and each consisted of 10 plants per treatment. The data was expressed as the mean ± SD of the three independent experiments.

The publicly available transcriptome data of TM-1 was used to assess the relative gene expression patterns (Zhang et al., 2015). RNA-Seq data were acquired and analyzed as described in the previous report (Yang et al., 2017; Qanmber et al., 2018). A heatmap was drawn by the MeV (Multiple Experiment Viewer) software (Saeed et al., 2003).

Four-day-old cotton seedlings were transferred into the MS medium containing DMSO and 5 μM AZD. The seedlings were collected at 24 hours post-treatment, frozen immediately in liquid nitrogen and stored at −80°C. Total RNA was extracted from samples using the RNAprep pure kit (Tiangen, Beijing) for RNA sequencing. Sequencing library preparation followed our previously described procedure (Song et al., 2017). Three biological replicates were sequenced using the Illumina Hiseq 2500 for each treatment. All raw sequencing data used in the current study have been deposited in the Sequence Read Archive database under accession number SRX2568809. The lncRNAs were identified from the novel and antisense transcripts of transcriptome assemblies by adopting the following steps: (1) the transcripts with low expression levels and with only a single exon were removed while the transcripts with two or more exons were selected; (2) the transcripts with length < 200 bp were removed; (3) the transcripts were searched against the Cuffcompare database to eliminate transcripts encoding proteins and protein coding domains; (4) the transcripts with FPKM (expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced) < 0.5 were removed; (5) the transcripts that failed to pass any one of the protein-coding-score tests using the Coding Potential Calculator (CPC), Coding-Non-Coding Index (CNCI) and Pfam Scan software were regarded as TUCP (transcripts of uncertain coding potential) (Kong et al., 2007; Sun et al., 2013; Finn et al., 2014). The optimized transcripts with no protein-coding potential were identified as lncRNAs and were listed in Supplementary Table 1.

The information about protein-coding genes and chromosomes was obtained according to the CottonGen database for cotton research (Yu et al., 2014). The program Circos was used to generate a diagram that showed the localization, abundance and expression level of lncRNAs and protein-coding genes in the cotton genome (Krzywinski et al., 2009).

The Cuffdiff software was used to identify the differentially expressed lncRNAs between DMSO control and AZD treatment (Trapnell et al., 2012). The fold changes were calculated via Log₂ (FPKM_AZD/FPKM_DMSO). LncRNAs with a |Log₂(Fold change)|≥1 and adjusted P ≤ 0.05 were selected as differentially expressed lncRNAs and displayed in Supplementary Table 2. The co-location protein-coding genes, which were upstream or downstream 100 kb away from lncRNAs, were used to analyze the GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) terms. The Pearson correlation coefficient was employed to explore the expression relationship between these lncRNAs and the neighboring protein-coding genes, and the lncRNAs with r_p > 0.95 were selected. GO analyses of the differentially expressed lncRNAs were performed by the GOseq program (Young et al., 2010). KOBAS software was used to test the enrichment of differentially expressed lncRNAs in KEGG pathways (Mao et al., 2005; Kanehisa et al., 2008).

Total RNA isolation and the generation of first-strand cDNA were performed as previously described (Song et al., 2017). The CFX96 system (BIO-RAD) was used to perform the QRT-PCR experiments by using an SYBR Premix Ex Taq II (TAKARA) kit. The primers used for QRT-PCR were designed by Primer Premier 5.0 and presented in Supplementary Table 3. Three biological replicates were performed and the reactions were performed in triplicate for each run. Expression levels of the target lncRNAs were normalized to GhHISTONE3 in cotton. The 2^−ΔΔCt analysis method was used to measure the relative expression levels of lncRNAs (Livak and Schmittgen, 2001). ANOVA analysis and a Tukey-Kramer multiple comparison test were used for data analysis. P < 0.05 was considered as statistically significant.

The TOR protein sequences of Arabidopsis were used as queries to search G. hirsutum database (Yu et al., 2014) using the blastp program with a cut-off E-value of 1e-5. Thereafter, four candidates were identified, and their sequences were extracted (Table 1). On account of the gene number and chromosome locus, the four GhTOR genes were designated as GhTOR1-At, GhTOR1-Dt, GhTOR2-At, and GhTOR2-Dt. Compared with the single TOR gene in diploid Arabidopsis, the number of TOR genes in allotetraploid cotton was significantly expanded. To dissect the evolutionary relationship of TOR genes between G. hirsutum and Arabidopsis, a phylogenetic tree was constructed by MEGA 5.0 software using the NJ method (Figure 1A). As shown in Figure 1A, GhTOR genes can be divided into two groups, and every group has two orthologs according to the phylogenetic tree, indicating that cotton TOR genes may have undergone duplication because only one AtTOR was found. Gene function is closely related with its structure; thus, we analyzed the gene structure of GhTORs (Figure 1B). Multiple-exons (50 to 66 exons) were displayed in GhTOR genes. We also found that GhTOR1-At and GhTOR1-Dt had highly similar exon patterns, and this was also the case between GhTOR2-At and GhTOR2-Dt. Similar domain organizations were also found between GhTORs and AtTOR proteins, except for GhTOR2-Dt, which lacked the FATC domain (Figure 1C). Collectively, these observations indicated that GhTORs have a close phylogenetic relationship with AtTOR and confirmed that TOR is also conserved in G. hirsutum.

AZD is a novel ATP-competitive specific inhibitor of the catalytic site of TOR kinase and has been successfully applied to research in the higher plants by independent groups (Montané and Menand, 2013; Dong et al., 2015; Li et al., 2015b; Xiong et al., 2017). AZD can directly bind to the kinase domain of TOR protein to inhibit the activity of both TORC1 and TORC2 (Chresta et al., 2010; Benjamin et al., 2011). AZD has shown relatively higher potency than other TOR inhibitors in both animals and plants (Montané and Menand, 2013; Dong et al., 2015). Thus, AZD was selected to dissect the biological functions of TOR in cotton (Figure 1D). Cotton seeds were germinated on the MS medium containing DMSO (as a control) or 5 μM AZD for 4, 6, and 8 days, respectively. These TOR-suppressed cotton seedlings showed a significant delay in the transition from heterotrophic to photoautotrophic growth compared with the cotton seedlings exposed to DMSO. The fresh weights and root lengths of TOR-suppressed cotton seedlings at 4, 6, and 8 days were only 86%, 60% and 69% of control seedlings, respectively (Figure 1E). Moreover, the number of lateral roots in the TOR-inhibited cotton seedlings also showed a significant decrease compared to that of the control seedlings. These data indicated that TOR plays a vital role in cotton growth.

To further investigate the cellular function of cotton TOR signaling, virus-induced gene silencing (VIGS) was used in this study (Liu et al., 2002; Pang et al., 2013). The conserved cDNA fragment of the four GhTOR genes was cloned into the virus vector pTRVRNA2 (pTRV2). The empty vector pTRV2 and vector pTRV2 containing the cDNA of GhPDS were used as negative and positive controls, respectively. These vectors were transformed into Agrobacterium cells and then inoculated into the cotyledons of the wild-type cotton CCRI24. As shown in Figures 2A–C, most of the growth of GhTOR-silenced plants was significantly retarded compared with that of the control seedlings. We found that most of the GhTOR-silenced plants exhibited pale yellow-green and shrinking leaves under normal growth conditions (Figure 2C), and the phenotype was consistent with the phenotype displayed in AZD-treated cotton seedlings (Figure 1D). A small part of GhTOR-silenced plants' leaves showed completely wilting phenotypes, which could be due to the lethal phenotype caused by the complete silence of GhTOR. The fresh weights of TOR-silenced cotton were also significantly reduced (Figure 2D). To validate that these phenotypes were caused by the reduced expression of the GhTOR genes, the expression levels of GhTOR genes were analyzed by QRT-PCR (Figure 2E). The top leaves of the cotton seedlings indicated in Figure 2A were harvested to perform QRT-PCR, respectively. Primers that targeted the conserved sequences of GhTOR1 (GhTOR1-At and GhTOR1-Dt) and GhTOR2 (GhTOR2-At and GhTOR2-Dt) were designed. The expression levels of GhTOR1 and GhTOR2 genes both decreased consistently in GhTOR-silenced plants. Our results confirmed that the retarded growth phenotype was caused by the reduced expression of the GhTOR genes.

Gene expression pattern is helpful to dissect gene biological function. To realize the biological function of TOR in cotton, we detected the expression patterns of the TOR genes using the publicly available transcriptome data of Gossypium hirsutum (Supplementary Figures 1A,B) (Li et al., 2015a). As shown in Supplementary Figure 1, the expression profiles of the orthologs between the GhTOR1-At, GhTOR1-Dt, and GhTOR2-Dt were similar in different tissues during seed germination, seed (fiber) development and root development, suggesting that their biological functions might be conserved. We found that all the GhTOR genes were widely expressed in the vegetative (root, stem and leaf) and reproductive (torus, petal, stamen, pistil, calycle and −3, −1, 0, 1, 3, 5, 10, 20, 25, and 35 days post-anthesis (DPA) ovule) tissues as well as in the fiber (5, 10, 20, and 25 DPA), suggesting that TOR genes have multiple biological functions in different tissues (Supplementary Figure 1A). GhTOR2_At/Dt have relatively low expression levels in cotton petals and stamens. GhTOR1_At/Dt were highly expressed in the stem, leaf and the late stage of ovule development (20 DPA and 25 DPA fibers). GhTOR1_At/Dt were highly expressed in seeds and roots at different sampling stages (Supplementary Figure 1B), indicating that GhTOR1_At/Dt play a potential role in the seed germination and root development.

LncRNAs can regulate gene expression through various mechanisms (Chekanova, 2015). The pair-end RNA-seq database has become a major resource to discover lncRNAs (Ilott and Ponting, 2013). To further understand TOR functions, paired-end RNA-Seq data of transcripts from CCRI24 and AZD-treated CCRI24 were used to find novel lncRNAs in Gossypium hirsutum seedlings. After data quality control, TOPHAT2 mapping and cufflinks transcript assembly, a total of 222,031 unique transcripts were assembled from the high-throughput RNA-Seq data. To distinguish lncRNA candidates, five sequential stringent filters to the transcripts were employed (Supplementary Figure 2). First, a total of 209,311 transcripts with two or more exons were selected, and 12,720 transcripts were removed. LncRNAs are longer than 200 nucleotides in length, so 1,840 transcripts that were shorter than 200 bp were removed, and 207,471 transcripts were recovered in the second step. The transcripts which were overlapped with the annotated exons in the Genecode database were removed by using the Cuffcompare software, and 23,106 transcripts were obtained. Next, the expression level of every transcript was quantified by the software Cuffquant, and a total of 17,064 transcripts with FPKM ≥ 0.5 were obtained. These transcripts were further filtered using CNCI, CPC, Pfam-scan and phyloCSF to predict the coding potential. Eventually, this pipeline provided 10,315 lncRNAs (Supplementary Figure 2A).

According to their genome locations, these lncRNAs were further classified into three types: long intergenic ncRNAs (lincRNAs), intronic ncRNAs (incRNAs), and natural antisense transcripts (NATs) (Supplementary Figure 2B). We found that in cotton, the majority of lncRNAs (9,356, 90.7%) were lincRNAs, and only a small portion (959, 9.3%) of the lncRNAs were lncNATs. Furthermore, no intronic lncRNAs were uncovered in our database. These observations are consistent with previous studies (Wang et al., 2015a). The previous studies in zebrafish, Arabidopsis, maize, Populus and tomato have shown that lncRNAs harbor fewer exons and are shorter than protein-coding genes (Pauli et al., 2012; Li et al., 2014a; Shuai et al., 2014; Wang et al., 2014a; Zhu et al., 2015). To further identify whether cotton lncRNAs share these features, the density of exon number, length and the open reading frame (ORF) length of cotton lncRNAs were analyzed compared with the cotton mRNA (Supplementary Figures 2C–E). We found that the majority of the lncRNAs only contained one to two exons. In contrast, the number of exons for the protein-coding genes ranged from one to twenty (Supplementary Figure 2C). There is no significant difference between the length of lncRNAs and mRNAs (Supplementary Figure 2D). Most of the ORF lengths of lncRNAs were < 250 nucleotides, while the majority of the ORF lengths of mRNA were >250 nucleotides (Supplementary Figure 2E). Taken together, compared with the protein-coding genes, the majority of the lncRNAs in cotton contain fewer exons, and are relatively shorter.

Next, we further analyzed the genomic landscape of cotton lncRNAs (Figure 3). A total of 10,315 lncRNAs were identified in the present study. 7,640 of these lncRNAs were mapped onto chromosomes of the recently released cotton reference genome (Yu et al., 2014). The remaining 2,670 lncRNAs were only mapped onto scaffolds. A Circos plot clearly showed that the cotton lncRNAs were widely transcribed from every chromosome. We found that lncRNAs and protein-coding genes shared similar density distributions, which are both lower in the pericentromeric heterochromatin regions than in the euchromatin. This result indicated that lncRNAs may share similar transcription patterns with the protein-coding genes, which is just the same as reported in the previous studies (Wang et al., 2015a). Furthermore, the expression levels of lncRNAs which mapped onto chromosomes were analyzed. LncRNAs with FPKM ≥0.5 were selected. A total of 528 lncRNAs were removed and 7,012 lncRNAs were obtained. Fold changes (FPKM_AZD/FPKM_DMSO) of these lncRNAs were calculated and are shown by the color of the lines on the inner track of Figure 3. The fold changes of 623 lncRNAs were ≥2 and are shown by red lines. The yellow lines represent 3,857 lncRNAs with fold change ≥1 and < 2. lncRNAs with fold change < 1 are depicted by green lines. The inner track provided a genomic landscape of lncRNA expression levels.

To identify the lncRNAs regulated by the TOR signaling pathway, the expression levels of lncRNAs were compared between DMSO control and AZD treatment in cotton seedlings. Following the procedure indicated in the methods, a total of 401 significantly differentially expressed lncRNAs were uncovered. Compared with the DMSO control, 296 of 401 lncRNAs were up-regulated in the cotton seedlings exposed to AZD, and the other 105 lncRNAs were downregulated (Figure 4A). To validate these differentially expressed lncRNAs, we randomly selected 10 of them, five were downregulated and five were upregulated. Their expression level changes indicated in RNA-seq data were validated by QRT–PCR experiments (Supplementary Figure 3). We also found that the expression patterns of these lncRNAs were consistent with that under TOR knockdown by VIGS.

Until now, the detailed mechanism of the interaction between mRNA and lncRNA is not clear. In previous studies, lncRNAs have been shown to preferentially locate next to genes with diverse functions (Dinger et al., 2008; Mercer et al., 2008; Guttman et al., 2009; Ponjavic et al., 2009; Cabili et al., 2011; Pauli et al., 2012; Wang et al., 2015a, 2018b). Thus, the protein-coding genes that are upstream or downstream 100 kb from the differentially expressed lncRNAs were selected to the Gene Ontology (GO) term analysis, and the top 10 enriched GO terms are shown in Figure 4B. The top two enriched GO terms were “defense response” and “response to biotic stimulus,” suggesting the potential function of the TOR signaling pathway in defense response through regulating cotton lncRNAs. The significantly differentially expressed lncRNAs potentially related to stress response are shown in Table 2, and a total of 21 lncRNAs were identified. We found that the lncRNAs that target the MLP-like protein 423 family, PYL and PP2C family may be regulated by TOR in cotton. KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis was also performed (Figures 4C,D). The top three enriched KEGG pathways of upregulated lncRNAs were “glutathione metabolism,” “arachidonic acid metabolism,” and “carotenoid biosynthesis.” The top three enriched KEGG pathways of down-regulated lncRNAs were “monoterpenoid biosynthesis,” “ribosome,” and “sphingolipid metabolism.”

As the phenotypes indicated in the TOR-inhibited and TOR-silenced cotton seedlings (Figures 1D, 2A), cotton seedling growths were significantly retarded. Plant cell growth is tightly linked to ribosome biogenesis and photosynthesis (Dong et al., 2015). The ribosome is a key component of cell growth control and protein synthesis. The “Ribosome” and “Ribosome biogenesis in eukaryotes” pathways were identified in the enriched downregulated KEGG pathways. LncRNAs involved in this process were listed in Supplementary Table 4A. A total of 37 downregulated lncRNAs were enriched in this pathway. Ribosomal protein (RP) is one of the important components of the ribosome (Ben-Shem et al., 2011). Twenty-nine of these downregulated lncRNAs targeted genes encoded ribosomal proteins, indicating the important role of TOR in regulating the cotton ribosome biogenesis. Transposons associate with the origins and transcribed regions of lncRNAs and eventually drive lncRNA expression (Ponting et al., 2009; Kapusta et al., 2013; Xu et al., 2015; Wang et al., 2016a,b; Sherafatian and Mowla, 2017). We analyzed the potential transposons that targeted these ribosome-related lncRNAs (Supplementary Table 4B). Transposons were identified in the sequences of LNC_000574, LNC_000642, LNC_004307, LNC_006480, LNC_006613 and LNC_007286 by WUBlastX. Most of the identified transposons were long terminal repeat retrotransposons (LTR), and three of them were LTR/copia transposons. These results indicated that transposon-related insertions may play a role in regulating lncRNA expression. One of the important functional patterns seen for lncRNAs is the interplay between miRNAs and lncRNAs, and lncRNAs could be targeted by miRNAs (Yoon et al., 2014). To identify whether lncRNAs are bona fide targets for miRNAs, we analyzed the selected lncRNAs associated with ribosomes using psRNATarget (Supplementary Table 4C) (Dai and Zhao, 2011). Twenty-four of the ribosome-related lncRNAs were targeted by miRNAs, and a total of 52 miRNAs were identified. In addition, one miRNA had multiple targets. For example, in the case of ghr-miR2949a-3p, it not only targeted LNC_000642 but also LNC_005100 and LNC_009692.

As one of the most important anabolic processes, photosynthesis also plays a vital role in plant cell growth. The TOR-inhibited and silenced cotton seedlings showed a significant delay in the transition from heterotrophic to photoautotrophic growth and exhibited pale yellow-green and shrinking leaves, which indicated the affected photosynthesis in cotton. The enriched KEGG pathways “carbon fixation in photosynthetic organism” and “photosynthesis” were also detected. Fifteen genes targeted by these down-regulated lncRNAs were found in these pathways and are displayed in Supplementary Table 5A. Genes encoded, including phosphoglycerate kinase (PGK1 and PGK3), alanine aminotransferase 2 (ALAAT1), aspartate aminotransferase 3 (ASP3), photosystem II reaction center protein K precursor, cytochrome c and photosynthetic electron transfer A, were found. Additionally, transposons and miRNAs were analyzed in the sequences of these lncRNAs. As displayed in Supplementary Table 5B, transposons were only identified in the sequences of LNC_004722. Furthermore, 26 miRNAs were found to target these photosynthesis-related lncRNAs, and eight of these lncRNAs were targeted by the miRNAs (Supplementary Table 5C). The above results further supported the phenotypes indicated in the TOR-inhibited cotton seedlings and the conserved functions of TOR as reported in other species (Martin et al., 2004; Ren et al., 2011; Dong et al., 2015).

As indicated in the GO terms analysis of differentially expressed lncRNAs, defense response-related GO terms were enriched. Environmental stresses significantly influenced plant growth and finally caused plant biomass losses, which were similar to the phenotypes indicated in the TOR-silenced cotton seedlings, including wilting and decreased fresh weights (Figure 2). Therefore, the expression patterns of GhTOR genes under abiotic stress conditions, i.e., cold, heat, salt, and drought were assessed, and the corresponding results were shown in Supplementary Figure 1C and Figure 5. The expression levels of GhTOR genes in Figure 5 were similar to the expression patterns shown in Supplementary Figure 1C. The expression levels of GhTOR genes were induced by diverse stresses, indicating that TOR responded to multiple stresses. When responding to different stresses, the expression levels of GhTOR1-At, GhTOR1-Dt, and GhTOR2-Dt increased; however, the transcript levels of GhTOR2-At were down-regulated. In response to cold, the expression levels of GhTOR1-At, GhTOR1-Dt, and GhTOR2-Dt reached a peak at 1 h, then maintained a relatively high and stable level. Meanwhile, in response to heat they reached a peak at 3 h. The expression levels of these three genes showed a gradual growth trend when responding to NaCl and PEG treatments. These results indicated that TOR genes may be involved and play potentially different roles in cotton stress response.

Expression pattern analysis of TOR genes indicated that TOR genes responded to different stress treatments (Figure 5 and Supplementary Figure 1). In our previous study, the differentially expressed protein coding genes of AZD-treated cotton seedlings were analyzed and identified (Song et al., 2017). The Pearson correlation coefficient was employed to explore the expression relationship between these lncRNAs and neighboring protein-coding genes, six (r_p > 0.95) of these lncRNAs were selected, including LNC_00512, LNC_00513, LNC_004695, LNC_004696, LNC_004697, and LNC_004698. LNC_00513 belonged to lincRNA, and the other five lncRNAs belonged to antisense lncRNAs. We found that the expression levels of the lncRNA-targeted genes Gh_A02G0500 and Gh_D02G0560, which both encoded MLP-like protein 423, were significantly upregulated in response to TOR inhibition, with the fold change of 14.8 and 5.1, respectively. As indicated in Table 2, Gh_A02G0500 was targeted by LNC_00512 and LNC_00513, and Gh_D02G0560 was targeted by LNC_004695, LNC_004696, LNC_004697, and LNC_004698. The expression levels of these above lncRNAs were downregulated under TOR inhibition with the fold changes varying from 0.2 to 0.4, indicating that TOR may be a positive regulator of these lncRNAs. Under TOR inhibition, these lncRNAs were downregulated, whereas their targets Gh_A02G0500 and Gh_D02G0560 were upregulated, suggesting the possibility that these lncRNAs may act as a negative regulator in the transcript process of these two genes. In the sequences of LNC_004698, six potential transposons were identified (Supplementary Table 6A) and all of them belonged to LINE/L1, indicating the potential role of transposons in regulating LNC_004698 expression. In the analysis of these stress related lncRNAs, we found 30 miRNAs targeted these lncRNAs, and 15 of these lncRNAs were targeted by miRNAs (Supplementary Table 6B). LNC_005731 and LNC_006174 were the target of ghr-miR399, which is engaged in salt tolerance of cotton (Deng et al., 2018).

To characterize the putative function of these lncRNA candidates in response to different stresses in Table 2, the expression levels of lncRNAs under the same stress treatments have been identified by QRT-PCR (Figure 6). We selected four lncRNAs that showed similar expression patterns with TOR genes for every stress treatment (Figure 6). In response to cold treatment, LNC_000512, LNC_004698, LNC_005100 showed a quick response and peaked at 1 h. In response to heat stress, LNC_004696, LNC_006174 and LNC_009095 reached the maximum at 3 h. We found that LNC_001977, LNC_004698, LNC_006174, and LNC_009095 displayed similar expression patterns when exposed to NaCl and PEG treatments, indicating that they may play a similar role in salinity and drought response. When responding to diverse stresses, all the transcript levels of LNC_007115 were downregulated, indicating that this lncRNA may be a negative regulator of stress response and may be regulated by GhTOR2-At, which showed a similar expression pattern. Gene expression profiles of these putative stress-related lncRNAs suggested that they may function as regulators downstream of TOR.

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