The susceptibility of sea-island cotton recombinant inbred lines to Fusarium oxysporum f. sp. vasinfectum infection is characterized by altered expression of long noncoding RNAs

Disease resistance is one of the most complicated yet important plant traits. The potential functions of long noncoding RNAs (lncRNAs) in response to pathogenic fungi remain unclear. In this study, we sequenced the transcriptomes of four different sea-island cotton (Gossypium barbadense) recombinant inbred lines (RILs) with susceptible, highly susceptible, highly resistant, or super highly resistant phenotypes and compared their responses to Fusarium oxysporum f. sp. vasinfectum (Fov) infection with those of their susceptible and resistant parents. Infection-induced protein coding genes were highly enriched in similar disease resistance-related pathways regardless of fungal susceptibility. In contrast, we found that the expression of a large number of Fov infection-induced lncRNAs was positively correlated with plant susceptibility. Bioinformatics analysis of potential target mRNAs of lncRNAs with both trans-acting and cis-acting mechanisms showed that mRNAs co-expressed or co-located with Fov-regulated lncRNAs were highly enriched in disease resistance-related pathways, including glutathione metabolism, glycolysis, plant hormone signal transduction, anthocyanin biosynthesis, and butanoate metabolism. Together these results suggest that lncRNAs could play a significant role in the response to pathogenic fungal infection and the establishment of disease resistance. The transcriptional regulation of these infection-susceptible lncRNAs could be coordinated with infection-susceptible mRNAs and integrated into a regulatory network to modulate plant-pathogen interactions and disease resistance. Fov-susceptible lncRNAs represent a novel class of molecular markers for breeding of Fov-resistant cotton cultivars.


Results symptoms of infected cotton.
Previous studies have shown that G. barbadense is more resistant to VW than G. hirsutum 20,35 but less resistant to Fov 36 . However, few studies have examined how G. barbadense responds to FW. Therefore, we obtained and crossed two G. barbadense cultivars with different levels of resistance to Fov infection, including a resistant cultivar (paternal plant, 06-146, [F]) and a susceptible cultivar (maternal plant, Xinhan-14, MX-14, [M]) 37 . In the F 6 RILs, four lines emerged based on their level of resistance, including two resistant (super highly resistant [SHR] and highly resistant [HR]) and two susceptible (highly susceptible [HS] and susceptible [S]) lines. In general, disease symptoms appear 4-7 days after inoculation. However, plant immune responses and transcriptional changes occur early during Fov infection 35,38 .
In order to profile plant transcriptional responses to Fov infection, seedlings at the one true-leaf stage were inoculated with Fov. Lower hypocotyls of cotton seedlings were then collected to assess responses to infection. A quantitative real-time PCR (qRT-PCR) assay was used to assess the expression of the Fov infection-related gene FOTG, and we found evidence of infection 28 h after inoculation ( Supplementary Fig. S1). We then extracted total RNA at a later time point, 40 h after inoculation, for RNA-seq experiments to assess the transcriptional response to Fov infection. We did not observe any symptoms 40 h after inoculation ( Supplementary Fig. S2). Previous studies have examined the growth of Fov in the lower hypocotyls of different plants within the first two days post infection (dpi). During this period of infection, green fluorescent protein-expressing strains of Fov were shown to grow on the root epidermis and adhere to the junctions of epidermal cells in tomato root 39 and melon 40 . In a study of Fov-infected cotton, it was shown that pathogens penetrate the roots 1 dpi, enter the vascular system, and move up through the plant to the hypocotyl vascular tissue 2-3 dpi, where infection causes vascular browning 38 . We observed vascular browning at a similar time point as this previous report. Based on these results, we anticipated that Fov infection would cause a transcriptional response in the lower hypocotyls 40 h after inoculation.  Table S1). Using TopHat2 alignment software 41 , 81.6% of all reads were successfully mapped to the current Upland cotton reference genome 34 . Global clustering analysis of Pearson's correlation coefficients between each sample pair revealed separate transcriptome patterns between control and infected groups (Fig. 1A). Differentially expressed gene (DEG) analysis (p value < 0.01 and |log 2 fold change (FC)| > 2) was performed to identify Fov-regulated genes in each pair of the infected and control RNA-seq samples for the two parental and four RIL lines. Clustering heatmaps of DEGs demonstrated a clear separation between control and infected groups, with most genes being downregulated (Fig. 1B). Upregulated and downregulated genes were consistent among infected and control genotypes, suggesting that the transcriptional response of mRNAs to Fov infection was highly conserved among the six genotypes (Fig. 1B). The number of DEGs was highest in the HS-FI versus HS-CK comparison (Fig. 1C). We detected 337 co-upregulated and 491 co-downregulated DEGs between the M-FI versus M-CK and F-FI versus F-CK comparisons, respectively, representing 26.79% and 22.15% of the upregulated and down-regulated DEGs, respectively ( Supplementary Fig. S3A). www.nature.com/scientificreports www.nature.com/scientificreports/ Analysis of the four RILs revealed similar results (Fig. 1D). These results suggested all the parental and RIL lines respond to Fov infection similarly, although substantial differences were presented as well.
EdgeR and other statistical pipelines were initially designed to assess RNA-seq data with replicates, although edgeR can also be used to assess RNA-seq data without replicates. Therefore, we used Chi-square test to analyze DEGs and compared the results with those of edgeR analysis. As shown in Supplementary Table S2, we found that the Chi-square test detected more differentially expressed (DE) mRNAs than edgeR in all groups except the HS group. The percentage of overlapping DEGs between edgeR and Chi-square tests was less than 30% of the total number of DEGs identified by the two methods, suggesting that differences exist between them.
We performed qRT-PCR experiments to validate the expression levels of DEGs. We selected six DEGs with elevated expression levels and one with reduced expression level in the HS line after infection. The Fov-regulated expression of all seven genes were well-detected by qRT-PCR analysis in the HS line ( Fig. 2A, Supplementary  Fig. S3B). We also found that Fov-upregulated expression of G1421, G1208, and G2549 were evident in almost all other groups. However, Fov-regulated expression of the other four protein-coding genes were either small or inversed in many other groups ( Fig. 2A, Supplementary Fig. S3B). www.nature.com/scientificreports www.nature.com/scientificreports/ To further explore the functions of DEGs, we performed functional clustering of upregulated genes and found similar enriched pathways among all four RILs ( Supplementary Fig. S4). Anthocyanin biosynthesis, phenylpropanoid biosynthesis, biosynthesis of secondary metabolites, and phenylalanine metabolism were highly enriched among the upregulated DEGs ( Supplementary Fig. S4). Similar Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were also enriched among downregulated DEGs, suggesting that expression of these genes was deregulated by Fov infection. DNA replication and plant hormone signal transduction were also enriched among downregulated DEGs ( Supplementary Fig. S4). We also obtained 452 co-upregulated DEGs in the two resistant groups and 473 co-upregulated DEGs in the two susceptible groups (Fig. 1D). Enriched KEGG pathways included glutathione metabolism, anthocyanin biosynthesis, flavonoid biosynthesis, biosynthesis of secondary metabolites, cysteine and methionine metabolism, and selenocompound metabolism (Fig. 2B). Many of these pathways have been associated with the response to fungal infection in plants. For example, flavonoid glycoside has been reported to exhibit antifungal activity against different Fusarium oxysporum f. sp. dianthi pathotypes 42 . Enriched BP terms in the Gene Ontology (GO) database also involved similar functions (Fig. 2C). Metabolic process, oxidation-reduction process, and transmembrane transport-related terms were most highly enriched (Fig. 2C), consistent with the important role of the oxidative burst in pathogen resistance.

Genome-wide identification and characterization of lncRNAs. A previous study revealed that lncR-
NAs play important roles in responding to VW in cotton 20 . Using an integrated approach, we identified 11,336 predicted lncRNAs from the 12 transcriptome datasets (Supplementary Table S3). Of these putative lncRNAs, 81.0% were long intergenic noncoding RNAs (lincRNAs), 14.4% were intronic lncRNAs, and 4.6% were natural antisense lncRNAs. In total, 63.9% and 36.1% of lncRNAs were distributed on the At and Dt subgenomes, respectively.
The global expression pattern of lncRNAs also differed between control and infected groups (Fig. 3A), demonstrating that Fov infection also greatly alters the expression pattern of lncRNAs. The mean length of lncRNA transcripts was shorter than that of protein-coding transcripts (1042.1 nt for lncRNAs and 3355.48 nt for protein-coding transcripts; Fig. 3B). The lengths of lncRNAs ranged from 201-89,370 nt, but more than 76.4% of www.nature.com/scientificreports www.nature.com/scientificreports/ lncRNAs were between 200 and 1000 nt in length (Fig. 3C). Consistent with previous studies 20, 21 , approximately 69.2% of lncRNAs had only one exon, and 30.8% had multiple exons (Fig. 3D). Calculated lncRNA expression levels were found to be lower than mRNA expression levels (t-test p value = 0, Fig. 3E). From these results, the predicted lncRNA profile was highly consistent with other studies, supporting the reliability of our prediction pipeline. From the boxplots of the 12 samples, no difference was found in the CK group, whereas a non-consistent expression pattern was found for the FI group (Fig. 3F). The two resistant samples showed downregulated expression after Fov infection (Fig. 3F). Inspection of global expression levels normalized to fragments per kilobase per million mapped reads (FPKM) for all mRNAs and lncRNAs indicated that the expression levels of most lncRNAs were lower than 10 FPKM (Fig. 3E,F).

exploration of pathogen-induced De lncRNAs in parental cultivars.
Because the two parental cultivars came from different genetic backgrounds and had different resistance levels, we hypothesized that they would demonstrate different responses to Fov infection. Using the edgeR 43 package of R software, 215 and 268 DE lncRNAs (false discovery rate [FDR] ≤ 0.05; |log2FC| ≥ 1) were found in the M-FI versus M-CK and F-FI versus F-CK comparisons, respectively. Although the number of DE lncRNAs was similar between these two groups, only 74 common DE lncRNAs were identified (Fig. 4A), indicating that lncRNA expression in response to Fov infection in these two cultivars was different. In both comparisons, Fov infection decreased expression of lncRNAs (Fig. 4B,C). By comparing DEGs and DE lncRNAs, we also detected more downregulated mRNAs in the M-FI versus M-CK and F-FI versus F-CK comparisons (Fig. 1C). As more DEGs and DE lncRNAs were detected in the paternal group, we concluded that this group was more sensitive to Fov infection.
The distribution of DE lncRNAs on chromosomes was also investigated. In M-FI versus M-CK and F-FI versus F-CK comparisons, DE lncRNAs were distributed on every chromosome. However, the lncRNA density (i.e., count by chromosome length) was not uniform across the two subgenomes. More DE lncRNAs were found on chromosomes A03, A05, A06, A12, D08, and D09 for the M-FI versus M-CK comparison (Fig. 4D), whereas more lncRNAs were found on chromosomes A01, A05, A11, A12, D09, and D11 for the F-FI versus F-CK comparison (Fig. 4D). This suggests that location bias plays a role in lncRNA induction by Fov infection.  www.nature.com/scientificreports www.nature.com/scientificreports/ respectively. In comparisons of the DE lncRNAs between these four groups, only 48 DE lncRNAs were found to be in common (Fig. 5A). These results suggest that different RILs respond to Fov infection in different ways. The number of lncRNAs regulated by Fov in the resistant strains HR (354) and SHR (444) was similar (Fig. 5B). However, the number of lncRNAs regulated by Fov in S (777) and HS (1768) RILs was much higher than in resistant RILs.

Identification of pathogen-induced DE lncRNAs in RILs.
In contrast to the DEGs shown in Fig. 1C, there were more upregulated than downregulated lncRNAs (Fig. 5B,C). The induced expression of a large number of Fov-responsive lncRNAs was positively correlated with susceptibility to Fov infection (Fig. 5B,C, boxed in blue). By performing Chi-square tests on expressed lncRNAs, www.nature.com/scientificreports www.nature.com/scientificreports/ we found fewer DE lncRNAs compared with edgeR. The percentage of overlapping DE lncRNAs varied from 13.25% to 39.62% (Supplementary Table S2).
We performed qRT-PCR experiments to verify the expression patterns of lncRNAs in the eight RIL samples. Thirteen upregulated DElncRNAs in response to infection from HS group were chosen for qRT-PCR experiments. The expression patterns of six selected genes are shown in Supplementary Fig. S5. The expression consistency of all 13 genes between their qRT-PCR and RNA-seq results was summarized in Supplementary Table S4. We observed 100% expression consistency for both Fov-susceptible RILs (HS and S), supporting the finding that lncRNAs are regulated by Fov infection. The expression consistency of these Fov-upregulated DElncRNAs was 92.3% and 46.2% consistency for the resistant RILs HR and SHR, respectively, suggesting the genotype-specific response of Fov infection.
We used k-means clustering to analyze DE lncRNAs and identified eight clusters. The k-means clustering results ( Supplementary Fig. S6 and Fig. 5D) corresponded with the results of hierarchical clustering shown in Fig. 6C. Sorting samples by decreasing Fov resistance levels, cluster 2, cluster 5, cluster 6, and cluster 8 revealed significantly elevated values after Fov infection (Fig. 5D and Supplementary Fig. S6). The reason that these four groups were not clustered together is that their expression values were not at the same level (Y axis in Fig. 5D), although their expression patterns were similar. Cluster 3, cluster 4 and cluster 7 contained 1483 DE lncRNAs and exhibited decreased expression in infected groups, corresponding with the downregulation of DE lncRNAs ( Supplementary Fig. S6A).
The distribution of DE lncRNAs on chromosomes was also investigated. In the four comparisons, Fov-induced DE lncRNAs were distributed across every chromosome (Fig. 5E). In the SHR-FI versus SHR-CK comparison, more DE lncRNAs were found on chromosomes A02, A03, A07, A08, A11, D11, and D12. In the HR-FI versus HR-CK comparison, more DE lncRNAs were found on chromosomes A03, A05, A11, A12, D02, and D05. In the S-FI versus S-CK comparison, more DE lncRNAs were found on chromosomes A05, A07, A11, D05, D08, and D11. In the HS-FI versus HS-CK comparison, more DE lncRNAs were found on chromosomes A05, A07, A11, D05, and D08. The non-uniform distribution of DE lncRNAs suggests that some genomic locations are enriched in DE lncRNAs. We performed clustering analysis to merge clusters of DE lncRNAs if the distance of adjacent lncRNAs was less than 10 kb. Circular plots revealed obvious clusters on chromosomes A05, A07, A11, D05, and D06 (Fig. 5F), which could serve as genomic hot spots for the response to fungal infection.

Functional enrichment analysis of pathogen-induced De lncRNAs. Regulated genes typically
show consistent or opposite expression patterns to their regulatory genes 44 and lncRNAs 45 . We performed co-expression analysis to predict the target mRNAs of DE lncRNAs. After filtering, 211,000 positive and 6607 negative pairs were retained (Fig. 6A). These filtered pairs included 4957 lncRNAs and 21,249 mRNAs, suggesting that potential lncRNA-mRNA regulatory networks are extensive. Of the DE lncRNAs shown in Fig. 5B, most DE lncRNA-mRNA pairs (96.03%) were positive. We then explored the features of each positive pair. The number of positive pairs per lncRNA was much higher for upregulated lncRNAs than for downregulated lncRNAs within the two susceptible RIL groups (paired t-test p value = 0.03, Fig. 6B), suggesting that some upregulated mRNAs may be regulated by upregulated lncRNAs in susceptible groups.
lncRNAs may function to modulate the transcription of nearby genes in a cis-acting manner 46 . Cis-acting regulation between lncRNAs and mRNAs also occurs. With a threshold distance of 10 kb between mRNAs and lncRNAs, 4394 lncRNA-mRNA cis-acting pairs were identified. Of these, 1173 pairs (26.70%) involving DE lncR-NAs were analyzed. Of these 1173 pairs, only 103 (8.78%) passed the co-expression threshold (Fig. 6C), and all were positively correlated, suggesting that there was no significant correlation in expression for cis-acting pairs and that lncRNAs may function in a trans-acting manner more often.
We then analyzed the functional pathways of the co-expressed and cis-acting mRNAs of DE lncRNAs. We focused on upregulated lncRNAs from the four clusters shown in Fig. 5D and Supplementary Fig. S6. First, we investigated the enriched KEGG pathways of mRNAs regulated by lncRNAs in a cis-acting manner. The most highly enriched pathway was the metabolic pathway (Fig. 6D). We also identified other enriched pathways in these four lncRNA clusters, including glutathione metabolism, plant hormone signal transduction, and SNARE interactions in vesicular transport (Fig. 6D). The same analysis was performed for mRNAs co-expressed with lncRNAs. The most highly enriched pathways were ribosome, glutathione metabolism, butanoate metabolism, and glycolysis/gluconeogenesis (Fig. 6E). Interestingly, glutathione metabolism appeared in both cis-acting and co-expression manners, suggesting that it plays an important role in lncRNA regulation in response to Fov infection.

Discussion
F. oxysporum is an important wilt pathogen that infects over 100 plant species, including cotton, in diverse ecological niches around the world. Breeding and cultivation strategies represent the most promising options to control this pathogen. Very few studies have been conducted to understand the mechanisms underlying Fov infection or to collect genetic marker information to facilitate molecular breeding of the important cotton variety sea-island cotton (G. barbadense). We previously created RILs from two parents selected from the field that showed divergent susceptibility to Fov infection. The RILs further diverged in their susceptibility, presenting an attractive opportunity to study sea-island cotton biomarkers of resistance and susceptibility to Fov.
In this study, we used a strand-specific RNA-seq approach to profile the transcriptome of parent and RIL plants with divergent susceptibility to Fov. We found that Fov-upregulated protein coding genes were highly enriched in known disease-resistance pathways, such as oxidation-reduction, glutathione metabolism, anthocyanin biosynthesis, and flavonoid biosynthesis pathways 6,13,[47][48][49][50] . However, downregulated genes were enriched in different biological functions, which indicates that susceptibility-determining protein coding genes were generally downregulated. Biological terms that were enriched in susceptibility-determining genes included ascorbate www.nature.com/scientificreports www.nature.com/scientificreports/ and aldarate metabolism. In summary, we obtained a large number of protein coding genes whose expression was responsive to Fov-infection and that could be used in molecular breeding in the future. Because of the lack RNA-seq replicates, the statistical power in identification of differentially expressed genes (DEGs) from RNA-seq www.nature.com/scientificreports www.nature.com/scientificreports/ data between the infected and control cotton for any individual group was limited. Nevertheless, we found that all the Fov-regulated genes showed highly similar expression response to Fov infection, and a large number of DEGs were overlapped, among multiple cotton genotypes. We therefore suggest the use of the overlapping DEGs as confident Fov-regulated genes for future study.
We used qRT-PCR to verify the Fov-regulated protein-coding (7) and non-coding genes (13) based on DEGs identified by edgeR analysis in HS group, demonstrating 100% consistency for both mRNA genes and lncRNAs in the same RIL genotype. The validation consistency was close to100% in some other RILs, but not in all. This finding agrees well with the presence of group variations in gene expression and in Fov response. Given the genetic difference among these six RILs, it is also possible that the discrepancy might reflect some RIL-specific differences in post-transcriptional regulation of these genes. For example, the presence of the cotton RIL-specific alternative splicing patterns could lead to altered ratio between splice isoforms which would not change the RNA-seq detected mRNA level, but might be sensitive to the specific qRT-PCR primers 51 .
We also performed Chi-square tests to analyze DEGs and DE lncRNAs and compared the results with those of edgeR. The low overlap in genes reflects the essential differences between these methods in calculating DEGs. Because the global expression level of lncRNAs was lower than that of mRNAs, we propose that the sensitivity of Chi-square tests is more dependent on gene expression level, as it detected fewer DE lncRNAs but more DEGs than edgeR analysis.
We also systematically identified sea-island cotton lncRNAs and studied their expression in response to Fov infection. In contrast to protein coding genes, lncRNAs showed extensive transcriptional responses to Fov infection in the two susceptible RIL plants, whereas resistant RIL plants showed similar transcriptional responses to the two parents. In the susceptible RILs, more lncRNAs were upregulated than downregulated and the number and magnitude of Fov-induced upregulated genes were both correlated with fungal susceptibility. Therefore, we identified a set of susceptibility-related lncRNAs. These findings provide additional molecular markers for breeding programs.
LncRNAs play important roles in various biotic and abiotic stress responses in plants by either regulating the transcriptional level of nearby genes in a cis-acting manner or by influencing other genes via trans-acting mechanisms 52,53 . Computational analysis of the two possible types of lncRNA regulation revealed interesting results. Firstly, co-expression analysis of potential trans-acting lncRNA regulatory mechanisms revealed that approximately 97% of lncRNA-mRNA pairs were positively co-expressed, which is much higher than in mammals 45 . It is possible that plant lncRNAs exert less negative and more positive regulation than mammalian lncRNAs. Secondly, protein coding genes co-expressed with DE lncRNAs are highly enriched in plant defense pathways, suggesting the importance of lncRNAs in the plant immune response. Thirdly, cis-acting regulation played a smaller role in modulating mRNA expression patterns, although the mRNAs adjacent to DE lncRNAs were enriched in plant defense-related pathways. Our results suggest that the two methods of lncRNA target prediction are of equal importance to explore lncRNA functions in plant defense. The predicted potential targets and functions of these lncRNAs will need to be verified in future studies.
Historically, a limited understanding of the genome-wide molecular responses of sea-island cotton to Fov infection and a lack of molecular markers related to Fov susceptibility have hindered the breeding of suitable high resistant cultivars. Transcriptome analysis of RIL plants with divergent Fov susceptibilities has revealed differences in mRNAs and lncRNAs that are correlated with fungal susceptibility. Moreover, lncRNA-mRNA network analysis revealed both positive and negative correlations between lncRNAs and their potential trans-acting and cis-acting mRNA partners that are related to disease resistance. Therefore, the findings of this study contribute to our understanding of disease resistance mechanisms in sea-island cotton and provide molecular markers to speed up breeding of disease-resistant cultivars.

Materials and Methods
plant material and pathogen inoculation. One highly aggressive strain of the defoliating fungus Fusarium oxysporum f. sp. vasinfectum from the College of Agronomy, Xinjiang Agricultural University was used for inoculation.
The seeds of a highly resistant G. barbadense cultivar (06-146) and a susceptible G. barbadense cultivar (Xinhan-14) were kindly provided by the Key Laboratory of Agricultural Biotechnology, Xinjiang Agricultural University. Four F 2:6 RILs were obtained and classified according to their level of resistance to Fov infection. These six cotton RILs were grown in germination boxes containing a sterilized mix of peat and sawdust. Before inoculation, we used sterilized scissor to cut fibril of root. Each seedling was inoculated with 10 mL of a suspension of 2 × 10 7 Fov spores per mL by watering injured roots of plants at the two-true-leaf growth stage. The inoculation was performed at room temperature (~25 °C), ≥60% humidity, 16 h/d light, and 8 h/d dark. After 40 h, susceptible cotton seedlings began to wilt. Control plants were not inoculated but were otherwise treated and sampled with distilled water in the same way. Hypocotyls were collected from seedlings for treatment after washing with 75% alcohol and sterile water. All samples were stored at −80 °C for future use.

RNA extraction and library construction. Total RNA was extracted using a Plant RNA EASYspin Plus
Kit (Aidlab, Peking, China) according to the manufacturer's instructions and treated with RQ1 DNase (Lot M6101, Promega, Madison, United States) to remove DNA. The quality and quantity of the purified RNA were determined by measuring the absorbance at 260 nm and 280 nm using a SmartSpec Plus (Bio-Rad, Hercules, CA, USA). RNA integrity was further verified by electrophoresis using a 1.5% agarose gel. To control for individual variation, we extracted total RNA from three randomly selected biological replicates for each sample and pooled the replicates together for further experiments.
www.nature.com/scientificreports www.nature.com/scientificreports/ A cDNA library was constructed using 10 μg of total RNA extracted from the hypocotyls of each G. barbadense sample. Polyadenylated RNA was purified and concentrated with oligo (dT)-conjugated magnetic beads (Invitrogen, Shanghai, China) before being used for directional RNA-seq library preparation. Purified mRNA was iron fragmented at 95 °C followed by end repair and 5′ adaptor ligation. Reverse transcription was then performed with RT primers harboring 3′ adaptor sequences and randomized hexamers. cDNA was purified and amplified, and 200-500 bp PCR products were purified, quantified, and stored at −80 °C until sequencing. cDNA clusters were generated and sequenced using an Illumina HiSeq2000 platform following the manufacturer's instructions to obtain 50-nt pair-end sequence reads. Two lanes were used for sequencing. Samples were separated by barcodes. qRt-pCR. In this study, to evaluate the validity of RNA-seq data, qRT-PCR was performed for the selected DEGs and DE lncRNAs. Polyadenylated RNA was enriched by oligo dT as in RNA-seq, which was then reverse transcribed into cDNA using M-MLV Reverse Transcriptase (Lot R011-01, Vazyme, Nanjing, China) and random primers. qRT-PCR was performed using the StepOne RealTime PCR System (QuantStudio ™ 6 Flex Real-Time PCR System Contains the OptiFlex ™ Optics System (Applied Biosystems ™ )) with the SYBR Green PCR Reagents Kit (Lot 11202ES08, Yeasen, Shanghai, China). The PCR conditions consisted of denaturing at 95 °C for 30 s, followed by 40 cycles of denaturing at 95 °C for 10 s and annealing and extension at 60 °C for 30 s. PCR amplifications were performed in triplicate for each sample and normalized using the cotton actin 7 gene as a reference gene. Cotton gene encoding Actin 7 (ACT7) was used as an internal control. Data were assessed using the comparative C T (ΔΔC T ) method 54 , and the relative expression level between the specific gene and ACT7 was thus obtained and presented. Primers for qRT-PCR analysis are listed in Supplementary Table S5. lncRNA identification and classification. By filtering low quality reads and adapters, high quality clean reads were obtained for downstream analyses. These reads were mapped to the Upland cotton genome 34 using TopHat2 41 software. Only mapped reads with one genomic location were used for analysis. FPKM was used to evaluate the expression level of genes and lncRNAs. Separated gene models from the same genotype were merged together using the Cuffmerge procedure 55 .
Novel transcripts were detected using Cufflinks and Cuffcompare 56 . Background noise was filtered from novel transcripts based on FPKM (>0.5), length (>200), coverage (>1), and status threshold (OK) 56 . The coding potential capability was calculated using the Coding Potential Calculator 57 (value < 0). Class code 'u' represented lincRNAs, 'x' represented long noncoding natural antisense transcripts (lncNATs), and 'i' represented intronic transcripts. lincRNA/protein-coding gene pairs were restricted to those that were non-overlapping and 1 kb away from protein-coding genes.
Identification of common and specific lncRNAs. All separate transcriptome gtf files of G. barbadense were merged into one gtf file using Cuffmerge with the parameter -g. These merged transcriptomes made it possible to compare the loci of lncRNAs from different sequencing samples using Cuffmerge. The class code 'u' represented specific lncRNAs among these sequencing samples. Furthermore, similar sequences were discarded to ensure the reliability of identified specific lncRNAs according to reciprocal BLASTN results with an E threshold value < 1e-10. The class code '=' represented core lncRNAs between G. barbadense that share fully equal loci. Reciprocal BLASTN analysis (E value < 1e-10) was also run to improve the confidence of identified core lncR-NAs, and only those with highly similar sequences were retained for further analysis.
lncRNA-mRNA co-expression analysis. Based on the expression of each mRNA and lncRNA, correlation coefficients and p values were obtained for each lncRNA-mRNA pair. We then filtered the results based on a threshold of absolute correlation coefficient no less than 0.95 and p value less than 0.01. In addition to positively correlated pairs, pairs with a correlation coefficient less than 0 were also included. These filtered gene pairs formed the expression network. expression and functional analysis. To measure the pattern of lncRNA expression and screen DE lncR-NAs, we applied the edgeR 43 package, which is specifically designed to analyze differential expression of genes using RNA-seq data. lncRNAs with FPKM < 0.5 in every sample were removed prior to analysis. To identify DE lncRNAs, FC (≥2 or ≤0.5) and FDR (≤0.05) cutoffs were applied.
To predict gene functions and calculate the functional category distribution frequency, KEGG analyses were performed using KOBAS bioinformatics resources 58 . Networks were constructed by calculating Pearson's correlation coefficients for DEGs. Cytoscape (v3.5.1) was used to display the co-expression network 59 . Circos software was used to plot the distribution of lncRNAs in the G. hirsutum genome 60 . statistical analysis. K-means was used to cluster the differently express pattern genes. Chi-square tests were also used to perform differential expression analysis. All values are presented as mean ± SD. Differences between means were determined using Student's t-tests. p values < 0.05 were considered statistically significant. All the statistical analysis was performed by R software.
Data deposition. Datasets supporting the results of this article are available in the NCBI Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) under Accession Number GSE95288.