Comparative physiology and transcriptome analysis allows for identification of lncRNAs imparting tolerance to drought stress in autotetraploid cassava

Background Polyploidization, pervasive among higher plant species, enhances adaptation to water deficit, but the physiological and molecular advantages need to be investigated widely. Long non-coding RNAs (lncRNAs) are involved in drought tolerance in various crops. Results Herein, we demonstrate that tetraploidy potentiates tolerance to drought stress in cassava (Manihot esculenta Crantz). Autotetraploidy reduces transpiration by lesser extent increasing of stomatal density, smaller stomatal aperture size, or greater stomatal closure, and reducing accumulation of H2O2 under drought stress. Transcriptome analysis of autotetraploid samples revealed down-regulation of genes involved in photosynthesis under drought stress, and less down-regulation of subtilisin-like proteases involved in increasing stomatal density. UDP-glucosyltransferases were increased more or reduced less in dehydrated leaves of autotetraploids compared with controls. Strand-specific RNA-seq data (validated by quantitative real time PCR) identified 2372 lncRNAs, and 86 autotetraploid-specific lncRNAs were differentially expressed in stressed leaves. The co-expressed network analysis indicated that LNC_001148 and LNC_000160 in autotetraploid dehydrated leaves regulated six genes encoding subtilisin-like protease above mentioned, thereby result in increasing the stomatal density to a lesser extent in autotetraploid cassava. Trans-regulatory network analysis suggested that autotetraploid-specific differentially expressed lncRNAs were associated with galactose metabolism, pentose phosphate pathway and brassinosteroid biosynthesis, etc. Conclusion Tetraploidy potentiates tolerance to drought stress in cassava, and LNC_001148 and LNC_000160 mediate drought tolerance by regulating stomatal density in autotetraploid cassava. Electronic supplementary material The online version of this article (10.1186/s12864-019-5895-7) contains supplementary material, which is available to authorized users.


Background
Cassava (Manihot esculenta Grantz) is a diploid plant (2n = 2× = 36) and an important cash and energy crop cultivated in Asia, Africa, and Latin America for its storage roots, making it critical for food security and economic development [1]. Water scarcity harms production when cassava is cultivated in severely waterdeficit regions, and although cassava can tolerant a wide range of adverse environmental conditions including drought, this can limit growth and survival [2].
Drought stress increases oxidative damage in plants [3] and reduces photosynthesis [4,5]. To cope with this, plants have evolved complex mechanisms such as deeper and thicker root systems [6], stomatal modulation systems including reduced aperture size and/or density to reduce water loss from transpiration [7,8], and accumulation of osmotic adjustment compounds [9]. At low to moderate concentrations, reactive oxygen species (ROS) such H 2 O 2 may act as second messengers in stress signalling. However, excessive H 2 O 2 production can trigger progressive oxidative damage to plant cells. Antioxidant enzymes scavenge excess H 2 O 2 and other ROS to protect plant cells from damage [10,11].
Doubling of the whole genome to generate polyploidy is ubiquitous among higher plant species, and the change to a polyploidy state increases abiotic stress tolerance in crop species. For example, tetraploidy in cenchrus, Arabidopsis and paulownia improves drought tolerance by lowering H 2 O 2 accumulation and enhancing ROS clearance [12][13][14]. Lower stomatal conductance and the abscisic acid (ABA) signalling pathway are involved in drought tolerance in the leaves of autotetraploid Rangpur lime (Citrus limonia) [15]. Autopolyploidization increases the potassium content and promotes tolerance to salinity stress in Arabidopsis [16]. Altered anatomy induced by polyploidy may confer stronger drought tolerance upon autotetraploid cassava [17]. However, the physiological and molecular advantages underlying these adaptations have not been widely investigated.
In the current study, we characterised physiological differences between diploid progenitor (2×) and autotetraploid (4×) cassava under standard and drought conditions. Comparative transcriptome analysis was used to identify differentially expressed genes (DEGs) and lncRNAs in the leaves of 2× and 4× cassava under wellwatered and dehydration conditions. Co-expression analysis revealed that two differentially expressed (DE) lncRNAs regulated six DEGs that improve drought tolerance in cassava by modulating stomatal density.

Results
Autotetraploid cassava displays stronger drought tolerance than diploid plants We exposed 'Xinxuan 048' 2× and 4× plants to soil with the same relative water content to compare their responses to drought stress. The RMSC was 30% for controls (Fig. 1a) and 4.5% for the drought stress treatment (Fig. 1b). After 15 days of withholding watering, at which time the RSMC was 4.5%, leaves of the 2× plants displayed moderate drooping, while leaves of 4× plants displayed only slight drooping (Fig. 1b). After withdrawing water for 30 days, followed by 2 days of recovery, all 2× plants were dead, while the growing points of 4× plants survived, and the upper leaves of autotetraploids remained green (Fig.  1c). Furthermore, the RWC of detached 4× leaves was higher than that of 2× plants (Fig. 2a). Analysis of water loss rate showed that 2× detached leaves lost water much more rapidly than those of 4× plants (Fig. 2b). Moreover, the transpiration rate of 4× seedlings was significantly lower than that of 2× seedlings (Fig. 2c). Thus, overall, 4× cassava plants were significantly more drought-tolerant than 2× plants.

Autotetraploidy alters drought-mediated stomatal function and photosynthetic capacity
Since water loss occurs mainly through stomatal movement, the reduced water loss in 4× plants prompted us to investigate the stomatal function of the two genotypes under control and 4.5% RSMC conditions. The stomatal density of 4× plants (7.24) showed a significant reduction compared with that of 2× plants (10.64) under control conditions, and the density of both genotypes increased with increasing drought stress, but the density of 4× plants (7.88) remained significantly lower than that of 2× plants (12.24) in dehydrated leaves (Fig. 3a). The average stomatal aperture in 2× plants under control conditions was 0.12 μm, compared with 0.09 μm in dehydrated 2× plants. The average stomatal aperture in control 4× plants was~0.16 μm, but this dropped 0.06 μm under drought condition (Fig. 3b). Thus, 4× leaves showed enhanced stomatal closure in response to drought, and stomatal conductance was markedly decreased in 4× leaves compared with 2× leaves under 4.5% RSMC (Fig. 3c). These results indicate that stress tolerance in 4× plants may be due to reduced transpiration rate via lesser extent increasing of stomatal density, smaller stomatal aperture size, and/or greater stomatal closure.
Next, we investigated photosynthesis parameters of the two genotypes. Photosynthesis is expected to be hindered due to lack of CO 2 under drought stress. Although the stomatal conductance of 2× plants was greater than that of 4× plants under drought stress, the net photosynthetic rate of 4× plants was significantly higher (Fig. 3d). This result is consistent with the higher SPAD value of 4× plants, although the difference was not significant (Fig. 3e).

Physiological effects of autotetraploidy
Six physiological traits were used to investigate dynamic changes in response to drought stress in 2× and 4× plants. An increase of more than 150% in H 2 O 2 production was observed in 2× leaves under drought stress, compared with less than 50% in drought-treated 4× leaves (Fig. 4a). By contrast, the relative increase in T-AOC and CAT activity following drought treatment was much less in 2× plants than 4× plants ( Fig. 4b and c). The MDA content is an indicator of the degree of lipid peroxidation, which can reflect damage to plant cell membranes [32]. Compared with controls, the MDA content was dramatically elevated in 2× dehydrated leaves but decreased in 4× dehydrated leaves, although the reduced value is not significant (Fig. 4d). The proline and soluble sugar content were also increased significantly in 2× plants with increasing drought stress, but were relatively unchanged in 4× plants ( Fig. 4e and f ). These results indicate that osmolyte biosynthesis is not influenced by tetraploidy. Overall, tetraploidy resulted in reducing accumulation of H 2 O 2 , and hence the ability to alleviate cell membrane injury, thereby increasing drought tolerance in 4× plants.

Identification of the lncRNAs and hierarchical clustering
In total, 1,232,088,086 bp raw reads were generated of 12 libraries by paired end sequencing. After moving the adapters and low-quality reads, 1,182,713,428 bp clean reads were mapped to the cassava genome. The mapped reads for each sample were assembled into transcripts using a reference-based approach. The correlation coefficient of the expression quantity of all transcriptome was more than 0.824 between each other among the 12 samples (Additional file 1: Figure S1). Expression levels of all coding genes and lncRNA transcripts were systematically estimated using FPKM values. In total, 2372 lncRNAs, including 821 antisense_lncRNAs and 1551 lincRNAs, were identified in this study (Additional file 2: Table S1). To calculate the degree of differential expression (DE) among lncRNAs, hierarchical clustering was performed using FPKM values of lncRNAs under control and drought stress conditions in 2× and 4× leaves. The results suggest that tetraploidization may have a limited effect on the mRNA transcriptome and lncRNA expression, since 2XCK and 4XCK clustered together, and 2XDR and 4XDR formed another cluster (Fig. 5a). Remarkably, significant DE was observed following drought   stress in both diploid and autotetraploid cassava, suggesting that DE lncRNAs may play an important role in responses to drought stress (Fig. 5b).

LncRNAs as potential miRNA precursors
By aligning miRNA precursors to the 2372 lncRNAs, we identified 18 lncRNAs as 21 known cassava miRNA precursors, including miR162, miR166c, miR408 and miR477c (Additional file 3: Table S2). A single lncRNA could serve as a precursor of several miRNAs, and different miRNAs could be targeted by the same lncRNA. Lnc_001314-miR171a was DE in response to drought.
Using a q-value < 0.05 and a fold change > 1 or < − 1 as thresholds, 1562 DEGs were found to be specifically responsive to drought in 4× sample (Fig. 6a), of which 687 genes were up-regulated and 875 were downregulated (Additional file 4: Table S3). A total of 5484 genes were commonly expressed ( Fig. 6a; Additional file 5: Table S4). We identified 2412 DEGs in 2× samples, including 1032 up-regulated and 1380 downregulated genes ( Fig. 6a; Additional file 6: Table S5). In 4× versus 2× plants subjected to dehydration, 814 DEGs were identified (Fig. 6b), including 288 down-regulated and 526 up-regulated DEGs (Additional file 7: Table S6). The reliability of the deep sequencing data was validated by quantitative real time PCR (qPCR) with gene-specific primers (Additional file 8: Table S7) for six randomly selected mRNAs (Additional file 9: Figure S2). The results showed that all the six genes displayed similar expression patterns in both RNA-seq and qPCR data. Notably, the number of drought-responsive mRNAs in 2× plants (2412) was higher than in 4× plants (1562), suggesting that 2× leaves might be more sensitive to drought than those of 4× plants, consistent with the more pronounced A closer inspection of the 5484 common genes identified a number of down-regulated genes. In particular, almost all the genes encoding subtilisin-like proteases, except MANES_05G206800 and MANES_06G139800, were more strongly down-regulated in 2× leaves than in 4× leaves following dehydration (Table 1). We used the qPCR to detect the expression of seven genes in the four samples, the results showed that MANES_05G206800 and MANES_06G139800 turned out to be opposite trend with the RNA-seq data, while the other five were consistent (Additional file 10: Figure S3). The discrepancy between qPCR and RNA-seq could be attributed to different statistical method. Consistently, a subtilisin-like protease gene mutant, sdd1-1, exhibited a two-to four-fold increase in stomatal density in Arabidopsis [33]. IiSDD1 was found to be autopolyploidy responsive and down-regulated by drought stress in autopolyploidy Isatis indigotica [34]. In the present work, the stomatal density was elevated in stressed 4× leaves (8.84%), but to a lesser extent than in 2× leaves (15.04%; Fig. 3a). Interestingly, the four subtilisinlike protease genes including, MANES_18G044300, MANES_18G039200, MANES_17G062600 and MANES_ 06G013200 were also found to be up-regulated expressed in 4× versus 2× plant subjected to dehydration (Additional file 7: Table S6). These results indicate a correlation between the regulation of stomatal density-related genes, resulting in a lesser extent of increasing stomatal density under dehydration stress, and consequently enhanced drought tolerance in 4× plants.
Among the 5484 common genes that were downregulated, genes involved in photosynthesis were particularly pronounced, including ferredoxin, ferredoxin-NADP reductase, beta-carbonic anhydrase, chlorophyll a-b binding protein, fructose-bisphosphate aldolase, photosynthetic NADH dehydrogenase subunit, and PsbP domaincontaining protein 3. Expression of these genes was reduced in dehydrated 4× leaves more than in 2× leaves ( Table 2), suggesting that photosynthesis was affected more in 4× plants than in 2× plants following drought treatment, consistent with the more pronounced reduction in stomatal aperture in 4× leaves under drought stress (Fig. 3b).
Regulating the expression of transcription factors (TFs) is an efficient strategy for amplifying transcriptional responses. As shown in Table 3, members belonging to the WRKY, MYB, and ERF TF families were represented. In particular, all seven ERF family members were downregulated in 4× plants under drought stress, and all 14 WRKY members except MANES_10G127100 and MANES_11G066500 were up-regulated under drought  (Table 3). Consistently, these three families of TFs appear to play important roles in drought stress signalling [35][36][37].
In order to explore the functions of the 86 lncRNAs specific to 4× plants, we constructed a co-expression   Tables S14 and S15). We used the qPCR to detect the expression of four lncRNAs with the corresponding trans target genes in the two samples, Me4XDR and Me4XCK. The results indicated that the expression is consistent with the RNAseq result (Fig. 7). Remarkably, we found that two of the DE lncRNAs might regulate six subtilisin-like protease DEGs; LNC_001148 (lincRNA) may regulate MANES_ 18G039200 (SBT1.7), MANES_06G139800 (SBT1.7), MANES_05G175300 (SBT1.6) and MANES_06G020400 (SBT2.5), while LNC_000160 (antisense lncRNA) may target MANES_16G094900 (SBT1.7) and MANES_ 18G044300 (SBT1.7). To validate the putative relationship between the two DE lncRNAs and six DEGs, expression levels were examined by qPCR. The results revealed similar lncRNA expression patterns to those obtained in the RNA-seq (Fig. 8), which suggests that the two lncRNAs identified using deep sequencing may target genes encoding subtilisin-like proteases via coexpression.
The functions of most of the 86 DE lncRNAs remain largely unknown, therefore, the genes co-expressed with 86 DE lncRNAs overlapping with DEGs were subjected to functional enrichment analysis using KEGG (Additional file 19: Table S16). The results indicated that 20 pathways were found responsive to drought stress. They included galactose metabolism, pentose phosphate pathway, plant hormone signal transduction, glycolysis/gluconeogenesis, biosynthesis of secondary metabolites and brassinosteroid biosynthesis, etc. (Additional file 20: Table S17).

Discussion
Autopolyploidy confers advantages over diploid counterparts, often manifested in enhanced adaptation to adverse environmental conditions. However, the molecular advantages remain largely unknown. Herein, we demonstrated the importance of autotetraploidy in response to drought in cassava. The 4× cassava plants potentiated stronger antioxidant and photosynthesis capacities than 2× plants, helping them cope better with drought stress. Autotetraploidy elevated the expression levels of two lncRNAs co-expressed with genes encoding subtilisin-like proteases regulating stomatal density under drought stress. Additionally, autotetraploidy decreased the stomatal density on the leaf epidermis, which in turn improves the RWC and reduces water loss, enhancing drought tolerance in cassava plants.
ROS participate in signal transduction, but are also toxic to cell membranes when present in excessive quantities, and can affect plant drought susceptibility [38]. We observed that H 2 O 2 in 4× stressed leaves were significantly lower than in 2× leaves (Fig. 4a). To detoxify excess drought-induced ROS, plants have developed a complex antioxidant system [39]. We found that 4× plants possessed a more efficient enzymatic antioxidant system involving CAT and T-AOC for reducing accumulation of H 2 O 2 under drought stress (Fig. 4b and c). MANES_05G1307001, homolog of catalase isozyme 1 in Ricinus communis, is up regulated in both of 2× and 4× stressed leaves compared with normal condition. However, the fold change of MANES_05G1307001 in 4× is more than that in 2×. Therefore, autotetraploidy may be associated with the regulation of antioxidation ability. Some antioxidant enzymes were found to be upregulated in stressed leaves, while others were downregulated, suggesting cellular redox status may be complex and dynamic. We found that levels of UDPglucosyltransferase were increased more or reduced less in leaves of 4× plants in response to drought compared with 2× leaves (Additional file 5: Table S4). Consistently, overexpression of Arabidopsis UGT79B2/B3 significantly enhances plant tolerance to drought stress by modulating anthocyanin accumulation, which enhances ROS scavenging [40]. Our results suggest that DEGs encoding UDP-glucosyltransferase may be tightly associated with improved drought tolerance in 4× by modulating ROS levels. Stomata are surrounded by two guard cells in the epidermis that regulate the shape and size of the pore aperture [41]. One of the earliest adaptive responses to drought in cassava leaves is stomatal closure and/or decreased stomatal density to reduce water loss [42]. Stomata control the balance between the uptake of CO 2 for photosynthesis and the release of water by modulating transpiration, thereby governing water use and abiotic stress tolerance [43]. It has been demonstrated that reducing the number of stomata does not affect carbon fixation due to increased CO 2 concentration [44]. Therefore, the decreased photosynthetic capacity of 4× plants is consistent with a lower transpiration rate, which results, at least in part, from enhanced stomatal closure during drought treatment (Fig. 3b-d). It is known that photosynthesis-related genes are down-regulated after drought treatment in many plants [45]. However, we found that the net photosynthesis rate of 4× plants remained higher than that of 2× plants under dehydration conditions, which might be attributed to higher SPAD values in 4× plants (Fig. 3e). Thus, 4× plants may maintain a higher photosynthetic rate, facilitating better adaptation to drought stress conditions by meeting increased energy demand.
Our transcriptomic data indicates that autotetraploidy influences the expression of genes encoding TFs involved in drought stress in cassava, and general response mechanisms integrating hormone signalling in response to external stimuli. TFs are efficient amplifiers of transcriptomic responses that regulate differential gene expression attributed to autotetraploidy. ERFs are key regulatory hub proteins in hormone and regulatory ROS-responsive gene expression that confer abiotic stress tolerance [37]. Some WRKY and MYB TFs are components of ABA-mediated stomatal movement involved in drought responses [36,46]. Three DE TF families were found to be in autotetraploid Arabidopsis under drought stress conditions [13]. A complex regulatory system controls drought tolerance in cassava. Most previous research in this area has Fig. 8 Comparison of expression results from RNA-Seq and qPCR methods for the two lncRNAs and co-expressed subtilisin-like protease family members. Data are means ± SD of three biological experiments. Cassava β-actin was used as an internal control focused on coding genes. LncRNAs are an important class of regulators in diverse biological processes involving complex mechanisms [47], and numerous lncRNAs have been identified in a few crop species including wheat (Triticum aestivum) [48], Medicago truncatula [49], Brassica napus [50], maize (Zea mays) [51] and cotton (Gossypium spp.) [52]. In our current study, we discovered 2372 lncRNAs, including 821 antisense_lncRNAs and 1551 lincRNAs, from 12 libraries. Li et al. [53] identified 682 lncRNAs from nine cassava samples. Differences between our current results and these previous results might be attributed to (i) differences in experimental design and transcriptome analysis, since the work of Li et al. [53] was based on 15-day-old seedling tissues under polyethylene glycol-simulated drought stress, and samples were collected from shoot apices and leaves; (ii) differences in bioinformatics strategies, since CPAT [54] and Pfam Scan [55,56] were employed to filter transcripts with coding potential and screen candidate lncRNAs. LncRNAs could execute their functions to respond to stress in either cis-acting or trans-acting in the genome via diverse mechanisms in plant [57,58]. Based on comparative transcriptome analysis, lncRNA16397 was found to be involved in Phytophthora infestans resistance by co-expression glutaredoxin in tomato (Lycopersicon esculentum Mill.) [59]. In the present study, LNC_001148 and LNC_000160, 898 bp and 2688 bp, respectively, were down-regulated by drought treatment in both 2× and 4× cassava plants. Co-expression analysis revealed that they appeared to regulate six genes encoding subtilisin-type proteinases. SDD1, belong to subtilisin serine proteinase family, which are known to negatively regulate stomatal density and distribution in Arabidopsis [33]. IiSDD1 participates not only in the drought responsive pathways, but also involves in autopolyploidy I. indigotica evolution [34]. Thus, relative to 2× dehydrated leaves, the higher expression of LNC_001148 and/or lower expression of LNC_ 000160 in 4× leaves may result in less down-regulation of target genes encoding subtilisin-type proteinases, and hence lesser extent of increasing stomatal density, thereby enabling 4× plants to better adapt to drought stress. Exactly how the two lncRNAs confer drought tolerance in cassava will be studied in future work.
In Table S14, LNC_000211 was down regulated with fold change more than 8 times following drought treatment in 4× plant specifically. LNC_000211 could trans regulate MANES_16G111000 (encoding Dehydrationresponsive element-binding protein), MANES_11G139300 (encoding ERF), MANES_18G039200 (encoding Subtilisin-like protease SBT1.7), MANES_08G148200 (encoding thioredoxin-like 1-1), respectively. Previous studies indicated that these four target genes were involved in drought tolerance in plants [33,37,60,61]. Therefore, LNC_000211 may play an important role in mediating the tolerance to drought in 4× cassava plant. Based on the KEGG enrichment analyses, we found that the trans target genes of the 86 DE lncRNAs were involved in galactose metabolism, pentose phosphate pathway, brassinosteroid biosynthesis, etc. The three pathways were reported to be associated with drought responsive in sugarcane (Saccharum officinarum L.), purging nut (Jatropha curcas) and Arabidopsis, respectively [62][63][64]. Taken together, our results suggest that 4× cassava implements divergent mechanisms to modulate the response to drought stress.
The current understanding of lncRNA regulation in response to drought stress is in its infancy in 4× cassava. These findings provide a comprehensive view of 4× DE lncRNAs, which will enable in-depth functional analysis.

Conclusion
Our study demonstrated that tetraploidy potentiates tolerance to drought stress in cassava. The co-expressed network analysis indicated that LNC_001148 and LNC_ 000160 in autotetraploid dehydrated leaves regulated six genes encoding subtilisin-like protease, thereby result in increasing the stomatal density to a lesser extent in autotetraploid cassava. This study helps to explain the role of autotetraploidy in conferring drought tolerance, and indicates the evolutionary potential of polyploidy in abiotic stress tolerance.

Water deficit treatment and water recovery
Cassava variety 'Xinxuan 048' (2× and 4×) used in this study were original from our lab [65]. Stem-propagated plants were grown in plastic pots (30 cm in height × 40 cm in diameter) containing well-mixed soil in March 2017. Each pot contained one cutting, and were placed in a greenhouse under a 16 h light/8 h dark photoperiod at the Guangxi Academy of Agricultural Sciences (GXAAS). Two-month-old cassava plants of each genotype were well-watered before drought stress treatment. Before dehydration treatment, each potted plant was watered until it was saturated to ensure consistency of water content. The moisture content of soil was measured by a AZS-100 soil moisture sensor (TRIME-PICO32, Germany). Control plants were well watered every 4 days. A relative soil moisture content (RSMC) of 30, 4.5% (after 15 days of withholding watering) and 1% (after 30 days of withholding watering until the top buds of diploid cassava displayed obvious wilting) was used for controls and drought stress treatments for two time points, respectively. After recovery for 2 days, the survival rate was determined and plants were photographed. All experiments were repeated in triplicate.

Measurement of water loss rate
Seventy-day-old plants were used to calculate the water loss rate, for which the fourth, fifth and sixth leaves (counting from the top of the plant) and petioles were excised from each plant. Three plants were included for each genotype. Detached leaves were placed on filter paper in a culture room under a 16 h light/8 h dark photoperiod. The abaxial surface was placed facing up for dehydration, and leaves were weighed every hour. Water loss was estimated from the percentage of fresh weight lost relative to the initial fresh weight. All experiments were repeated in triplicate.

Relative water content (RWC) assay
Four plants each of each genotype were used to measure the RWC. Seventy-day-old plants were detached at the fifth leaf (counting from the top of the plant) and weighed immediately (M1), then placed in water for 18 h under dark conditions, dried using a filter paper, then weighed again (M2). Saturated leaves were finally ovendried for 24 h at 65°C to a constant weight (M3), and the RWC was measured using eq. 1: Drought stress treatment and plant sampling One hundred-day-old cassava plants were used for drought stress analysis. Before dehydration treatment, each potted plant was watered until it was saturated to ensure consistency of water content. The moisture content of soil was measured by a AZS-100 soil moisture sensor (TRIME-PICO32, Germany). Control plants were well watered every 4 days. RSMC of 30 and 4.5% (after 15 days of withholding watering) was used for controls and drought stress treatments, respectively, and all treatments included three plants per genotype. The fifth leaf of each plant was used to measure stomatal aperture size, photosynthetic capacity, and physiological indices. The lower epidermis of leaves was used for examination as described previously by Zhou et al. [65]. The length and width of stomata were measured using 200 guard cells, and the number of stomatal openings was counted in 10 microscopic fields using three independent replicates. A LI-6400 portable photosynthesis system was used to measure the net photosynthetic rate, transpiration rate, and stomatal conductance. The relative chlorophyll content in the fifth leaves of each plant was determined by a soil plant analysis development (SPAD) value measured by an SPAD-502 Plus chlorophyll meter (Konica Minolta, Japan). Physiological indices, namely catalase (CAT) activity, total antioxidant capacity (T-AOC), H 2 O 2 content, and malondialdehyde (MDA) content, were estimated using reagent kit cat. # A007-2, A015, A064 and A003 (Institute of Nan Jing Jian Cheng Bioengineering), respectively. Proline content and total soluble protein content were measured using the method of Bates et al. [66] and Guy et al. [67], respectively.
Plants used for RNA-seq were treated as described above, and parallel leaves were frozen in liquid nitrogen and stored at − 80°C.

Identification of DEGs and DE lncRNAs
Cuffdiff (v2.1.1) was used to calculate fragments per kilobase of transcript per million mapped reads (FPKM) values of both lncRNAs and coding genes in each sample [70]. Gene FPKMs were computed by summing the FPKMs of transcripts in each gene group. To assess the three biological replicates, the log 2 (fold_change)-transformed FPKM values were performed. Cuffdiff provides statistical routines for determining DE in digital transcript using a model based on the negative binomial distribution [70]. A default q-value < 0.05 was set as the threshold for DE. Genes with log 2 (fold_change) values > 0 were deemed up-regulated, while genes with log 2 (fold_change) values < 0 were considered down-regulated.

Trans target genes prediction of lncRNAs
Trans targets genes of lncRNAs were identified from expression correlations between lncRNAs and coding genes using custom scripts. The pearson correlation coefficient method was used to calculate the correlation between lncRNA and mRNA among samples. Correlations corresponding to a coefficient > 0.95 or < − 0.95 were selected for the functional enrichment analysis.