Wild rice, including Oryza nivara and Oryza rufipogon, which are considered as the ancestors of Asian cultivated rice (Oryza sativa), possess high genetic diversity and serve as a crucial resource for breeding novel cultivars of cultivated rice. Although rice domestication related traits, such as seed shattering and plant architecture, have been intensively studied at the phenotypic and genomic levels, further investigation is needed to understand the molecular basis of phenotypic differences between cultivated and wild rice. Drought stress is one of the most severe abiotic stresses affecting rice growth and production. Adaptation to drought stress involves a cascade of genes and regulatory factors that form complex networks. O. nivara inhabits swampy areas with a seasonally dry climate, which is an ideal material to discover drought tolerance alleles. Long noncoding natural antisense transcripts (lncNATs), a class of long noncoding RNAs (lncRNAs), regulate the corresponding sense transcripts and play an important role in plant growth and development. However, the contribution of lncNATs to drought stress response in wild rice remains largely unknown.
Here, we conducted strand-specific RNA sequencing (ssRNA-seq) analysis of Nipponbare (O. sativa) and two O. nivara accessions (BJ89 and BJ278) to determine the role of lncNATs in drought stress response in wild rice. A total of 1246 lncRNAs were identified, including 1091 coding–noncoding NAT pairs, of which 50 were expressed only in Nipponbare, and 77 were expressed only in BJ89 and/or BJ278. Of the 1091 coding–noncoding NAT pairs, 240 were differentially expressed between control and drought stress conditions. Among these 240 NAT pairs, 12 were detected only in Nipponbare, and 187 were detected uniquely in O. nivara. Furthermore, 10 of the 240 coding–noncoding NAT pairs were correlated with genes enriched in stress responsive GO terms; among these, nine pairs were uniquely found in O. nivara, and one pair was shared between O. nivara and Nipponbare.
We identified lncNATs associated with drought stress response in cultivated rice and O. nivara. These results will improve our understanding of the function of lncNATs in drought tolerance and accelerate rice breeding.
Rice is one of the most important crops in the world and a major source of food for billions of people. Asian cultivated rice (Oryza sativa) was domesticated from wild rice species, including Oryza rufipogon and Oryza nivara, and experienced a bottleneck effect that severely reduced its genetic diversity  and decreased its viability in the natural environment . In contrast to cultivated rice, wild rice species possess higher genetic diversity, which contributes to its greater resistance to biotic and abiotic stresses, and this characteristic of wild rice is crucial for understanding and improving the stress tolerance of cultivated rice.
Drought stress is one of the most severe abiotic stresses affecting crop yield. In the natural environment, plants adapt to drought stress by employing various strategies, such as speeding up their life cycle to avoid drought stress, reducing water loss, and improving water use efficacy . Drought tolerance is a complex trait involving the regulation of a number of physiological and biochemical processes, including stomatal density , leaf rolling , osmotic adjustment , and root system development , at different development stages. These mechanisms of drought stress response involve genes belonging to various families including WRKY, MYB, NAC, ABRE, PP2C, and SnRK2 [8,9,10,11,12,13].
Long noncoding RNAs (lncRNAs) constitute a large fraction of the transcriptome that does not encode proteins , and play important roles in various biological processes, such as genome stability , vernalization , telomere maintenance , transcriptional activation , and other developmental processes . Various transcriptome studies have revealed that pervasive transcription from noncoding transcripts can give rise to functional lncRNAs . lncRNAs are classified as long intergenic noncoding RNAs (lincRNAs) and long noncoding natural antisense transcripts (lncNATs) . lncNATs are transcribed from the opposite strand of sense RNA in the same genomic regions, and may regulate the expression of sense RNA [22,23,24]. For example, lncNATs affect genes on the opposite strand and increase starch content and grain weight in rice . Furthermore, lncNATs have been studied at the genomic level in many species, including human (Homo sapiens) , mouse (Mus musculus) [27, 28], rice [29,30,31,32], Arabidopsis thaliana [21, 33], maize (Zea mays) , Plasmodium falciparum , and yeast (Saccharomyces cerevisiae) . In particular, lncNATs are correlated with the response to various abiotic and biotic stresses [29, 34, 37]. For example, a recent study in Arabidopsis revealed that NATs affect plant thermotolerance .
Although our understanding of drought stress response in plants has improved substantially, this topic needs to be investigated further, given the ongoing global climate change. Firstly, more genetic variability that contributes to drought tolerance needs to be identified in the plant germplasm. Secondly, since most of the previous studies focused on cultivated rice, the drought stress response in wild rice remains largely unknown and needs special attention.
Here, we performed strand-specific RNA sequencing (ssRNA-seq) analysis of Nipponbare (O. sativa ssp. japonica; hereafter referred to as Nip) and two accessions (BJ89 and BJ278) of O. nivara collected from its native habitats of Cambodia and Laos. O. nivara inhabits swampy areas with a seasonally dry climate . We identified 1091 coding–noncoding NAT pairs, of which 240 pairs were differentially expressed between control and drought stress conditions, and 187 pairs were specifically found in O. nivara. Furthermore, according to the GO enrichment analysis of sense transcripts, 10 coding–noncoding NAT pairs were correlated with drought stress, of which nine pairs were uniquely identified in BJ89 and BJ278. Thus, we identified numerous lncNATs correlated with drought stress in O. nivara. These lncNATs potentially play important roles in the response to drought stress, and will provide new insights into the mechanism of drought tolerance in O. nivara, thus facilitating breeding the cultivated rice varieties.
ssRNA-seq and transcript assembly
To examine lncNAT expression patterns in wild rice under drought stress treatment, one cultivated rice accession (Nip) and two O. nivara accessions (BJ89 and BJ278) were grown under control and drought conditions (Fig. 1a–c). Leaves of BJ89 and BJ278 showed lower water loss than those of Nip (Fig. 1d). Furthermore, the two O. nivara accessions exhibited a higher survival rate than Nip after 25 days of drought stress treatment (Fig. 1e). These physiological data suggest that O. nivara accessions BJ89 and BJ278 are more drought tolerant than Nip at the seedling stage.
Next, we conducted ssRNA-seq analysis of these three accessions treated with or without 25% polyethylene glycol (PEG-4000; w/v) for 10 days. A total of 18 strand-specific cDNA libraries were constructed from leaf tissues, with three replicates per accession in both control and drought stress treatments. In total, 445.5 million paired-end reads (2 × 125 bp) were generated by ssRNA-seq using Illumina HiSeq 2500, of which 373.6 million reads (83.9%) mapped perfectly on to the Nip reference genome (Table S1). Pearson correlation coefficients of the three biological replicates of each accession were greater than 0.9, indicating the high reproducibility of our ssRNA-seq data (Fig. S1). Among the 373.6 million paired-end reads, we identified 62,201 transcripts, including 17,583 novel transcripts (8905 known gene loci and 2704 new gene loci) with hisat2  and stringtie  (Fig. 2a). We also conducted principal component analysis (PCA) of gene expression data. The results showed that PC2 clearly distinguished between the control and drought treated samples, while PC3 separated the different accessions (Fig. S2).
Identification of lncRNAs
To identify lncRNAs, novel transcripts larger than 200 nt were mapped against the Rfam 13.0 database to exclude micro RNAs (miRNAs), ribosomal RNAs (rRNAs), and other small noncoding RNAs . Then, any transcripts with a coding potential, according to Coding Potential Calculator (CPC)  and Pfam with HAMMER scan , were filtered out (Fig. 2a). Finally, a total of 1246 lncRNAs were identified, including 940 in Nip, 959 in BJ89, and 974 in BJ278 (reads in the same accession under both control and drought stress conditions were combined to identify lncRNAs) (Fig. 2b and Table S2). Among the 1246 lncRNAs, 692 (55.5%) were common to all three accessions, 306 (24.6%) were uniquely found in at least one of the two O. nivara accessions, and 111 (8.9%) were present in both O. nivara accessions (Fig. 2b). The expression profiles of lncRNAs were more different between the three accessions than between the drought and control conditions of the same accession (Fig. 2c).
Of the 1246 lncRNAs, a total of 394 lncRNAs were differentially expressed between control and drought stress treatments (118 in Nip, 227 in BJ89, and 174 in BJ278). Among these 394 lncRNAs, 23 were common to all three accessions; 34 were identified as being shared between one of the O. nivara accessions and Nip; 45 were shared between the two O. nivara accessions (BJ89 and BJ278); and 139, 92, and 61 were uniquely found in BJ89, BJ278, and Nip, respectively (Fig. 3a, b).
Identification of lncNATs and NAT pairs
Based on their location relative to the gene coding regions, 675 of the 1246 lncRNAs were long intergenic noncoding RNAs (lincRNAs; lncRNAs located in intergenic regions), and 571 lncRNAs were long noncoding natural antisense transcripts (lncNATs; lncRNAs overlapped with coding genes on the opposite DNA strand) (Table S2). The lincRNAs contained fewer exons than lncNATs; 72.7% lincRNAs contained only one exon (Fig. S3). By contrast, mRNAs contained more exons than both lincRNAs and lncNATs (Fig. S3). In addition, mRNAs showed higher expression variation than lincRNAs and lncNATs under either control or drought stress condition at the genome level (Fig. S4).
It has been shown that lncNATs can regulate the expression of sense transcripts [45, 46], and each strand of a NAT pair could potentially represent a protein-coding gene. Therefore, in addition to identifying NAT pairs from lncRNAs, we scanned NAT pairs across the whole genome of the three accessions, based on the annotation of the Nip reference genome. All transcripts annotated in the Nip reference genome were integrated and assembled with the ssRNA-seq data generated in this study. A total of 8529 NAT pairs with overlapping regions greater than 25 nt were identified according to a previous study . According to the coding capacity of the sense–antisense pair , 86.88% (7410) of NAT pairs were coding–coding pairs (both transcripts with protein-coding capacity), 0.33% (28) were noncoding–noncoding pairs (both transcripts represented lncRNAs), and 12.79% (1091) were coding–noncoding pairs (one strand showed protein-coding capacity, while the other strand represented an lncRNA) (Table S3), and each transcript (transcripts with protein-coding capacity or lncRNAs) could flank a few different transcripts. Depending on the direction and location of the sense and antisense transcripts, 61.33% of the 8529 NAT pairs were enclosed (one transcript fully embedded in the other), 23.38% were divergent (head-to-head, 5′-end overlap), and 15.29% were convergent (tail-to-tail, 3′-end overlap) (Fig. 3d).
Of the 8529 NAT pairs, 5866 were detected through ssRNA-seq analysis of the three accessions, of which 4783, 4813, and 4596 were detected in Nip, BJ89, and BJ278, respectively (Fig. 3c, Fig. S5, Table S3). Of the 5866 NAT pairs, 62.2% (3651) were shared among all three accessions (2792 coding–coding pairs, 13 noncoding–noncoding pairs, and 846 coding–noncoding pairs) (Fig. S5, Table S3); 8.8% (517) were expressed only in Nip (462 coding–coding pairs, 5 noncoding–noncoding pairs, and 50 coding–noncoding pairs) (Fig. S5, Table S3); and 18.5% (1083) were uniquely expressed in O. nivara accessions (998 coding–coding pairs, 8 noncoding–noncoding pairs, and 77 coding–noncoding pairs) (Fig. S5, Table S3).
NAT pairs responsive to drought stress
To clarify the response of NAT pairs to drought stress, we first determined the differences in gene expression patterns between control and drought treatments. In detail, we identified differentially expressed genes (DEGs) using the following criteria: fold change (FC) ≥ 2.0 and false discovery rate (FDR) ≤ 0.01. A total of 3934 (4110 transcripts), 5880 (6235 transcripts), and 5036 (5294 transcripts) DEGs were identified between control and drought treatments in Nip, BJ89, and BJ278, respectively (Fig. 3a, Table S4). To identify genes within biological processes related to drought stress, GO enrichment analysis was performed on all DEGs (FC ≥ 2.0 and FDR ≤ 0.01) and highly differentially expressed genes (HDEGs) (FC ≥ 4.0 and FDR ≤ 0.01) identified in each accession. Based on all DEGs, a total of 57 GO terms in the three accessions, and most terms were related to primary metabolic pathways essential for plant growth and development, such as ‘biosynthetic process’, ‘cellular biosynthetic process’, and ‘primary metabolic process’ (Fig. S6). In addition, different GO terms were enriched in the three accessions in response to drought; for example, the ‘response to water’ GO term was uniquely enriched in BJ278 (Fig. S6).
Based on the HDEGs (FC ≥ 4.0 and FDR ≤ 0.01), 63 GO terms in the biological process category were enriched, including 10 terms related to stress, such as ‘response to jasmonic acid stimulus’, ‘oxidation reduction’, and ‘gibberellin metabolic process’ [47, 48] (Fig. 4a, Table S5). Among these ten stress related terms, three (‘response to chemical stimulus’, ‘response to stimulus’, and ‘response to stress’) were detected in both BJ89 and Nip; three terms (‘jasmonic acid mediated signaling pathway’, ‘response to jasmonic acid stimulus’, and ‘response to biotic stimulus’) were uniquely enriched in BJ89; and four terms (‘response to abiotic stimulus’, ‘oxidation reduction’, and ‘response to water and gibberellin metabolic process’) were only enriched in Nip. In BJ278, only one GO term (‘carbohydrate metabolic process’) was enriched. These results suggest that Nip, BJ89, and BJ278 employ different mechanisms to respond to drought stress (Fig. 4a, Fig. S6). A total of 134 HDEGs were enriched in these 10 stress related GO terms. Among these 134 genes, 48 were found only in O. nivara (either one or both accessions); 12 were found only in Nip; 35 were shared between Nip and one of the two O. nivara accessions (Table S5). In addition, there are 39 genes that are differentially expressed in all three accessions between control and drought treatment. However, these 39 genes exhibit more stronger response in O. nivara (Table S5). Antisense transcription could silence or concordantly regulate the sense transcripts . To detect NAT pairs responsive to drought stress, both sense and antisense transcripts of each NAT pair showing differential expression (FC ≥ 2.0 and FDR ≤ 0.01) between control and drought stress conditions were identified as differentially expressed NAT pairs. A total of 369 differentially expressed NAT pairs were identified, of which 240 were coding–noncoding pairs (193 in BJ89, 96 in BJ278, and 53 in Nip). Additionally, among these 240 differentially expressed NAT pairs, 23 were common in all accessions; 18 were shared between Nip and BJ89 or BJ278; 12 were present only in Nip; and 187 were found only in O. nivara accessions (Fig. S7, Table S6).
According to the effect of antisense transcripts on sense transcripts, we classified the differentially expressed NAT pairs into two categories, as described previously : discordant (sense and antisense transcripts showing opposite expression patterns) and concordant (sense and antisense transcripts expressed coordinately). A total of 24 discordant NAT pairs were identified (one shared among all accessions; two shared between BJ89 and BJ278; one shared between Nip and BJ278; 6, 6, and 8 uniquely found in Nip, BJ89, and BJ278, respectively) (Fig. S8), including one NAT pair discordant in Nip but concordant in BJ89, and one pair discordant in BJ89 and BJ278 but concordant in Nip (Fig. 4b and Table S6). Additionally, 102 concordant NAT pairs were upregulated (30 in Nip including one downregulated in BJ89, 61 in BJ89, and 54 in BJ278) (Fig. 4b). Among the upregulated concordant NAT pairs, 10 were shared among all accessions, 7 were shared between Nip and BJ278, 13 were found only in Nip, and 72 were found only in O. nivara (Fig. S9 and Table S6). A total of 245 concordant NAT pairs were downregulated (45 in Nip, 205 in BJ89, and 82 in BJ278), of which 17 were shared among all three accessions, 14 were shared between Nip and BJ89 or BJ278, 14 were found only in Nip, and 200 were found only in O. nivara (Fig. 4b, Fig. S10, and Table S6). RT-qPCR validated that the ssRNA-seq results based on the two samples (Nip and BJ278) for two NAT pairs that could be validated based on RT-qPCR using sequence specific primers (Fig. S11).
Among the 10 coding–noncoding NAT pairs related to GO enrichment terms for the drought stress response, nine were uniquely found in O. nivara, and one was common to all accessions (Table 1 Fig. S12, see Fig. S13 for the sequence alignment of the 10 NAT pairs). Furthermore, among the nine coding–noncoding NAT pairs uniquely found in O. nivara, six were correlated with response to stress, three were correlated with the jasmonic acid stimulus pathway, and one was correlated with oxidation reduction (Table 1).
Antisense transcripts are present in various organisms and play important roles in regulating gene expression. For example, the gene encoding the famous transcriptional repressor FLOWERING LOCUS C (FLC), which delays flowering time in Arabidopsis, is repressed at warm temperatures by COLD INDUCED LONG ANTISENSE INTRAGENIC RNA (COOLAIR), an antisense RNA, via histone demethylation [46, 49, 50]. Antisense lncRNAs can also upregulate gene expression; for example, an lncRNA transcribed from promoter region of the Pcdhα gene leads to DNA demethylation of the CTCF binding sites and the activation of sense promoters . Natural cis-antisense transcripts also define the function of short interfering RNAs (siRNAs) and affect their biogenesis; for example P5CDH and SRO5 regulate salt tolerance by generating two types of siRNAs in Arabidopsis . In addition, NATs also contribute to heterochromatin formation and DNA methylation, and suppress gene expression in tumorous cells . However, the function of conserved lncRNAs can vary across different species. For example, in human and mouse, the conserved lncRNAs exhibit different subcellular localization patterns .
O. nivara, one of the wild progenitor species of cultivated rice, inhabits swampy areas with a seasonally dry climate. O. nivara possesses greater genetic diversity than cultivated rice and thus is a valuable resource for breeding rice cultivars with desirable traits, such as stress resistance. In this study, we identified 240 coding–noncoding differentially expressed NAT pairs, of which both of the coding genes on one strand and the lncRNAs on the opposite strand were differentially expressed between control and drought stress conditions. Furthermore, we identified 187 coding–noncoding NAT pairs in the leaves of O. nivara accessions, of which 10 were correlated with drought stress. These findings are expected to facilitate further analysis of drought stress tolerance in wild rice. However, to identify more drought resistance genes and NAT pairs, a more comprehensive analysis of different tissues and developmental stages of wild rice should be conducted. It is interesting that even the two accessions from O. nivara have different NAT pairs. The variation could be resulted from the genetic variation of within this species, which has a much high genetic variation than cultivated O. sativa.
Furthermore, these NAT pairs need to be validated through functional analysis, and the functions of a vast majority of NAT pairs remain unknown. Therefore, further studies are needed to understand the molecular mechanism of action of these coding–noncoding NAT pairs and their regulatory network under drought stress conditions.
In this study, we performed ssRNA-seq analysis of cultivated rice Nipponbare and two O. nivara accessions, and systematically identified numerous genes and NAT pairs responsive to drought stress in all three accessions. Overall, the results of this study will enhance our understanding of the evolution of NATs and drought tolerance in rice, thus facilitating rice breeding.
Plant material and growth conditions
Oryza nivara accessions BJ89 and BJ278 were collected from Cambodia and Laos, respectively, and were deposited at the PE Herbarium, Institute of Botany, Chinese Academy of Sciences. All experiments were conducted in an environmentally-controlled growth chamber (Percival Scientific, Inc.). Seeds were incubated at 42 °C for at least 5 days to break dormancy. Seeds were then placed on a wet filter paper in a culture dish at 30 °C for 3 days. The most uniformly geminated seeds were sown in a bottom-less 96-well plate placed in a dish containing Kimura B culture solution , and incubated in the growth chamber at 30 °C day/25 °C night temperature under a 12-h light/12-h dark photoperiod. The culture solution was renewed every 3 days. To conduct drought stress and control treatments, 12-day-old seedlings were transferred to a Kimura B culture solution containing 25% PEG-4000 (w/v) or no PEG-4000, respectively, and incubated for 10 days. Three biological replications were performed for each accession in each treatment.
Strand-specific cDNA library preparation and ssRNA-seq
Total RNA was extracted from leaves of five plants that had been subjected to 10-day-drought or control treatments using Promega Total RNA Isolation System (Z3100). Then, mRNA was fragmented in the fragmentation buffer (Ambion). The first strand of cDNA was reversed transcribed with SuperScript II (Invitrogen) using random primer (TAKARA). Next, dTTP was replaced by dUTP for the synthesis of the second strand. The repaired double-stranded cDNA was ligated with Illumina TruSeq adaptor and digested using USER enzyme (NEB). Finally, 350–450 bp fragments were recovered and purified, and 18 cDNA libraries (3 accessions × 2 treatments × 3 replicates) were sequenced on the Illumina HiSeq 2500 platform.
Measurement of water loss and survival rate
A water loss experiment was conducted to investigate the variation in the rate of water loss rate from leaves among Nip, BJ89 and BJ278. Five topmost fully expanded leaves were collected from 25-day-old seedlings of each accession and weighted every 30 min at room temperature. Three replications were performed for each accession. The relative water loss rate was calculated using the following equation:
where FW and CW represent the fresh and current weight of leaves, respectively.
To examine the differences in survival rate among the three accessions, 12-day-old seedlings of each accession were treated with Kimura B culture solution containing 25% PEG-4000 (w/v) for 4 weeks and then transferred to Kimura B culture solution without PEG for 1 week. The survival rate of each accession was calculated as the ratio of alive plants to all plants.
Transcript assembly and lncRNA identification
Paired-end (2 × 125 bp) strand-specific reads of each accession were combined and mapped to the Nip reference genome sequence (IRGSP-1.0) using the hisat2 (2.1.0) software . Stringtie was used to assemble transcripts for regions with read coverage greater than 5X . Then, transcripts of the three accessions were integrated using Cuffmerge and Cuffcompare utilities in the Cufflinks package . Gene annotations of the Nip reference genome (http://rapdb.dna.affrc.go.jp) were integrated, and gene expression levels (estimated as FPKM [Fragments Per Kilobase of transcript per Million mapped reads] values) and DEGs were determined using Cuffdiff in the Cufflinks package . Transcripts longer than 200 nt were mapped against the Rfam 13.0 database to exclude miRNAs, rRNAs, and other small noncoding RNAs . Then, Coding Potential Calculator (CPC)  and Pfam  were used to filter out potential coding transcripts. Gene function annotations were available from the RAP-DB (http://rapdb.dna.affrc.go.jp/index.html). GO enrichment analysis was conducted using agriGO .
RNA was extracted from 14-day-old seedlings using the SV Total RNA Isolation System (Promega). First-strand cDNA was obtained using RevertAid First Strand cDNA Synthesis Kit (ThermoFisher). RT-qPCR (one cycle of 95 °C for 15 s, followed by 45 cycles of 95 °C for 15 s, 56 °C for 45 s) was performed using the qTOWER3 Real-Time PCR Thermal Cycler (Analytik Jena, Germany) with TB Green Premix Ex TaqII (TaKaRa, Japan) according to the manufacturers’ instructions. The 2-ΔΔCT method was used to determine the gene expression level. The expression level of all transcript was normalized against actin. Each experiment was repeated with three independent biological replicates, and RT-qPCR reactions were performed with three biological replicates and three technical replicates for each sample. The primers used in this experiment are listed in Table S7.
The ssRNA-seq datasets generated during the current study are available in the China National Center for Bioinformation (https://bigd.big.ac.cn/databases) under the accession number: CRA003736. The samples used in this study were deposited at the PE Herbarium, Institute of Botany, Chinese Academy of Sciences.
Ishimaru K, Shirota K, Higa M, Kawamitsu Y. Identification of quantitative trait loci for adaxial and abaxial stomatal frequencies in Oryza sativa. Plant Physiol Bioch. 2001;39(2):173–7. https://doi.org/10.1016/S0981-9428(00)01232-8.
Yue B, Xue W, Xiong L, Yu X, Luo L, Cui K, et al. Genetic basis of drought resistance at reproductive stage in rice: separation of drought tolerance from drought avoidance. Genetics. 2006;172(2):1213–28. https://doi.org/10.1534/genetics.105.045062.
Barbez E, Dunser K, Gaidora A, Lendl T, Busch W. Auxin steers root cell expansion via apoplastic pH regulation in Arabidopsis thaliana. Proc Natl Acad Sci U S A. 2017;114(24):E4884–93. https://doi.org/10.1073/pnas.1613499114.
Miao J, Li X, Li X, Tan W, You A, Wu S, et al. OsPP2C09, a negative regulatory factor in abscisic acid signalling, plays an essential role in balancing plant growth and drought tolerance in rice. New Phytol. 2020;227(5):1417–33. https://doi.org/10.1111/nph.16670.
Mao H, Wang H, Liu S, Li Z, Yang X, Yan J, et al. A transposable element in a NAC gene is associated with drought tolerance in maize seedlings. Nat Commun. 2015;6(1):8326. https://doi.org/10.1038/ncomms9326.
Tao Z, Kou Y, Liu H, Li X, Xiao J, Wang S. OsWRKY45 alleles play different roles in abscisic acid signalling and salt stress tolerance but similar roles in drought and cold tolerance in rice. J Exp Bot. 2011;62(14):4863–74. https://doi.org/10.1093/jxb/err144.
Yoshida T, Fujita Y, Sayama H, Kidokoro S, Maruyama K, Mizoi J, et al. AREB1, AREB2, and ABF3 are master transcription factors that cooperatively regulate ABRE-dependent ABA signaling involved in drought stress tolerance and require ABA for full activation. Plant J. 2010;61(4):672–85. https://doi.org/10.1111/j.1365-313X.2009.04092.x.
Gong ZZ, Xiong LM, Shi HZ, Yang SH, Herrera-Estrella LR, Xu GH, et al. Plant abiotic stress response and nutrient use efficiency. Sci China Life Sci. 2020;63(5):635–74. https://doi.org/10.1007/s11427-020-1683-x.
Hu WL, Jin L, Xu A, Wang YF, Thorne RF, Zhang XD, et al. GUARDIN is a p53-responsive long non-coding RNA that is essential for genomic stability. Nat Cell Biol. 2018;20(4):492–502. https://doi.org/10.1038/s41556-018-0066-7.
Feretzaki M, Pospisilova M, Valador Fernandes R, Lunardi T, Krejci L, Lingner J. RAD51-dependent recruitment of TERRA lncRNA to telomeres through R-loops. Nature. 2020;587(7833):303–8. https://doi.org/10.1038/s41586-020-2815-6.
Grossi E, Raimondi I, Goni E, Gonzalez J, Marchese FP, Chapaprieta V, et al. A lncRNA-SWI/SNF complex crosstalk controls transcriptional activation at specific promoter regions. Nat Commun. 2020;11(1):936. https://doi.org/10.1038/s41467-020-14623-3.
Coudert AE, Pibouin L, Vi-Fane B, Thomas BL, Macdougall M, Choudhury A, et al. Expression and regulation of the Msx1 natural antisense transcript during development. Nucleic Acids Res. 2005;33(16):5208–18. https://doi.org/10.1093/nar/gki831.
Wang H, Chung PJ, Liu J, Jang IC, Kean MJ, Xu J, et al. Genome-wide identification of long noncoding natural antisense transcripts and their responses to light in Arabidopsis. Genome Res. 2014;24(3):444–53. https://doi.org/10.1101/gr.165555.113.
Wang KC, Yang YW, Liu B, Sanyal A, Corces-Zimmerman R, Chen Y, et al. A long noncoding RNA maintains active chromatin to coordinate homeotic gene expression. Nature. 2011;472(7341):120–4. https://doi.org/10.1038/nature09819.
Cabianca DS, Casa V, Bodega B, Xynos A, Ginelli E, Tanaka Y, et al. A long ncRNA links copy number variation to a polycomb/trithorax epigenetic switch in FSHD muscular dystrophy. Cell. 2012;149(4):819–31. https://doi.org/10.1016/j.cell.2012.03.035.
Katayama S, Tomaru Y, Kasukawa T, Waki K, Nakanishi M, Nakamura M, et al. Antisense transcription in the mammalian transcriptome. Science. 2005;309(5740):1564–6. https://doi.org/10.1126/science.1112009.
Oono Y, Yazawa T, Kanamori H, Sasaki H, Mori S, Matsumoto T. Genome-wide analysis of rice cis-natural antisense transcription under cadmium exposure using strand-specific RNA-Seq. BMC Genomics. 2017;18(1):761. https://doi.org/10.1186/s12864-017-4108-5.
Zhou X, Sunkar R, Jin H, Zhu JK, Zhang W. Genome-wide identification and analysis of small RNAs originated from natural antisense transcripts in Oryza sativa. Genome Res. 2009;19(1):70–8. https://doi.org/10.1101/gr.084806.108.
Lu T, Zhu C, Lu G, Guo Y, Zhou Y, Zhang Z, et al. Strand-specific RNA-seq reveals widespread occurrence of novel cis-natural antisense transcripts in rice. BMC Genomics. 2012;13(1):721. https://doi.org/10.1186/1471-2164-13-721.
Zhang YC, Liao JY, Li ZY, Yu Y, Zhang JP, Li QF, et al. Genome-wide screening and functional analysis identify a large number of long noncoding RNAs involved in the sexual reproduction of rice. Genome Biol. 2014;15(12):512. https://doi.org/10.1186/s13059-014-0512-1.
Liu J, Jung C, Xu J, Wang H, Deng S, Bernad L, et al. Genome-wide analysis uncovers regulation of long intergenic noncoding RNAs in Arabidopsis. Plant Cell. 2012;24(11):4333–45. https://doi.org/10.1105/tpc.112.102855.
Xu J, Wang Q, Freeling M, Zhang X, Xu Y, Mao Y, et al. Natural antisense transcripts are significantly involved in regulation of drought stress in maize. Nucleic Acids Res. 2017;45(9):5126–41. https://doi.org/10.1093/nar/gkx085.
Siegel TN, Hon CC, Zhang Q, Lopez-Rubio JJ, Scheidig-Benatar C, Martins RM, et al. Strand-specific RNA-Seq reveals widespread and developmentally regulated transcription of natural antisense transcripts in Plasmodium falciparum. BMC Genomics. 2014;15(1):150. https://doi.org/10.1186/1471-2164-15-150.
Yassour M, Pfiffner J, Levin JZ, Adiconis X, Gnirke A, Nusbaum C, et al. Strand-specific RNA sequencing reveals extensive regulated long antisense transcripts that are conserved across yeast species. Genome Biol. 2010;11(8):R87. https://doi.org/10.1186/gb-2010-11-8-r87.
Borsani O, Zhu J, Verslues PE, Sunkar R, Zhu JK. Endogenous siRNAs derived from a pair of natural cis-antisense transcripts regulate salt tolerance in Arabidopsis. Cell. 2005;123(7):1279–91. https://doi.org/10.1016/j.cell.2005.11.035.
Li Y, Li X, Yang J, He Y. Natural antisense transcripts of MIR398 genes suppress microR398 processing and attenuate plant thermotolerance. Nat Commun. 2020;11(1):5351. https://doi.org/10.1038/s41467-020-19186-x.
Cai Z, Zhou L, Ren NN, Xu X, Liu R, Huang L, et al. Parallel speciation of wild rice associated with habitat shifts. Mol Biol Evol. 2019;36(5):875–89. https://doi.org/10.1093/molbev/msz029.
Kalvari I, Argasinska J, Quinones-Olvera N, Nawrocki EP, Rivas E, Eddy SR, et al. Rfam 13.0: shifting to a genome-centric resource for non-coding RNA families. Nucleic Acids Res. 2018;46(D1):D335–42. https://doi.org/10.1093/nar/gkx1038.
Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35(Web Server issue):W345–9.
Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, Mitchell AL, et al. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 2016;44(D1):D279–85. https://doi.org/10.1093/nar/gkv1344.
Sun Q, Csorba T, Skourti-Stathaki K, Proudfoot NJ, Dean C. R-loop stabilization represses antisense transcription at the Arabidopsis FLC locus. Science. 2013;340(6132):619–21. https://doi.org/10.1126/science.1234848.
Marquardt S, Raitskin O, Wu Z, Liu F, Sun Q, Dean C. Functional consequences of splicing of the antisense transcript COOLAIR on FLC transcription. Mol Cell. 2014;54(1):156–65. https://doi.org/10.1016/j.molcel.2014.03.026.
Liu F, Quesada V, Crevillen P, Baurle I, Swiezewski S, Dean C. The Arabidopsis RNA-binding protein FCA requires a lysine-specific demethylase 1 homolog to downregulate FLC. Mol Cell. 2007;28(3):398–407. https://doi.org/10.1016/j.molcel.2007.10.018.
Liu F, Marquardt S, Lister C, Swiezewski S, Dean C. Targeted 3′ processing of antisense transcripts triggers Arabidopsis FLC chromatin silencing. Science. 2010;327(5961):94–7. https://doi.org/10.1126/science.1180278.
Canzio D, Nwakeze CL, Horta A, Rajkumar SM, Coffey EL, Duffy EE, et al. Antisense lncRNA transcription mediates DNA demethylation to drive stochastic protocadherin α promoter choice. Cell. 2019;177(3):639–53 e615. https://doi.org/10.1016/j.cell.2019.03.008.
Yu W, Gius D, Onyango P, Muldoon-Jacobs K, Karp J, Feinberg AP, et al. Epigenetic silencing of tumour suppressor gene p15 by its antisense RNA. Nature. 2008;451(7175):202–6. https://doi.org/10.1038/nature06468.
Guo CJ, Ma XK, Xing YH, Zheng CC, Xu YF, Shan L, et al. Distinct processing of lncRNAs contributes to non-conserved functions in stem cells. Cell. 2020;181(3):621–36 e622. https://doi.org/10.1016/j.cell.2020.03.006.
Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5. https://doi.org/10.1038/nbt.1621.
We thank members of the Guo lab for providing helpful suggestions for this study.
This work was supported by the National Natural Science Foundation of China (31925004 to Y.-L.G.), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB27010305 and XDA08020103 to Y.-L.G.), The Innovative Academy of Seed Design, Chinese Academy of Sciences (Y.-L.G.), and China Postdoctoral Science Foundation (2020 M680749 to Y.-C.X.).
Yong-Chao Xu and Jie Zhang contributed equally to this work.
Authors and Affiliations
State Key Laboratory of Systematic and Evolutionary Botany, Institute of Botany, Chinese Academy of Sciences, Beijing, 100093, China
Yong-Chao Xu, Jie Zhang, Dong-Yan Zhang, Ying-Hui Nan, Song Ge & Ya-Long Guo
University of Chinese Academy of Sciences, Beijing, 100049, China
Yong-Chao Xu, Jie Zhang, Dong-Yan Zhang, Ying-Hui Nan, Song Ge & Ya-Long Guo
Y.-L.G. conceived the study; Y.-C.X., J.Z., S.G., and Y.-L.G. analyzed the data; J.Z. and Y.-H.N. performed phenotypic measurements; D.-Y.Z. collected the leaf samples and extracted RNA for RNA-seq; S.G. provided the O. nivara accessions; Y.-C.X. and Y.-L.G. wrote the manuscript, with contributions from all authors. The author(s) read and approved the final manuscript.
With the permission to collect, two accessions of Oryza nivara accessions were identified and collected by Professor Song Ge from Cambodia and Laos, respectively, and now deposited at the PE Herbarium (BJ89 and BJ278), Institute of Botany, Chinese Academy of Sciences. The study comply with relevant institutional, national, and international guidelines and legislation.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Pearson correlation coefficient of 18 strand-specific RNA-seq (ssRNA-seq) datasets. Gene expression levels (estimated as FPKM [Fragments Per Kilobase of transcript per Million mapped reads] values) were used for this analysis. Three replicates were performed for each accession. NC and NP indicate Nipponbare (cultivated rice; Nip) samples under control and drought stress conditions, respectively; BJ278C and BJ278P represent BJ278 (O. nivara; wild rice) samples under control and drought stress conditions, respectively; BJ89C and BJ89P represent BJ89 (O. nivara; wild rice) samples under control and drought stress conditions, respectively. Fig. S2. Principal component analysis (PCA) of FPKM data of Nip, BJ89, and BJ278 obtained under drought stress and control treatments. NC and NP indicate Nip samples under control and drought stress conditions, respectively; BJ278C and BJ278P represent BJ278 samples under control and drought stress conditions, respectively; BJ89C and BJ89P represent BJ89 samples under control and drought stress conditions, respectively. Fig. S3. Number of exons in long intergenic noncoding RNAs (lincRNAs), long noncoding NATs (lncNATs), and mRNAs. Fig. S4. Expression levels (log2FPKM) of different RNAs in Nip, BJ89, and BJ278. NC and NP indicate Nip samples under control and drought stress conditions, respectively; BJ278C and BJ278P represent BJ278 samples under control and drought stress conditions, respectively; BJ89C and BJ89P represent BJ89 samples under control and drought stress conditions, respectively. Red triangles represent coefficient of variation. Fig. S5. Venn diagram of NAT pairs expressed in Nip, BJ89 and BJ278. Fig. S6. Gene Ontology (GO) enrichment analysis based on differentially expressed genes (DEGs) under drought stress. Most GO terms were related to primary metabolic pathways. Colors indicate P-values of GO terms. Fig. S7. Venn diagram of differentially expressed NAT pairs under drought stress in Nip, BJ89 and BJ278. Fig. S8. Venn diagram of discordant NAT pairs under drought stress in Nip, BJ89, and BJ278. Fig. S9. Venn diagram of upregulated concordant NAT pairs under drought stress in Nip, BJ89, and BJ278. Fig. S10. Venn diagram of downregulated concordant NAT pairs under drought stress in Nip, BJ89, and BJ278. Fig. S11. The expression level of randomly selected two concordant NAT pairs that could design strand-specific primers. The relative expression level of Os02t0258800–01 (a, sense transcript) and MSTRG.6860.1 (b, antisense transcript) were both increased after PEG treatment. The relative expression level of BJ278 were increased in both Os02t0504000–01(c, sense transcript) and MSTRG.7513.1 (d, antisense transcript) after treatment but the transcripts in Nip showed no difference between control and treatment. “***” indicated that the p-value was less than 0.001. Fig. S12. Venn diagram of coding–noncoding NAT pairs of which enriched in GO terms related to drought stress in Nip, BJ89, and BJ278. Fig. S13. Sequence alignments of 10 coding–noncoding NAT pairs related to GO enrichment terms for the drought stress response. The top sequence was sense transcript and the sequence below was antisense transcript. Table S1. Statistics of paired-end reads generated by ssRNA-seq analysis of 18 cDNA libraries of Nipponbare (Nip), BJ89, and BJ278 samples treated with or without drought stress. Table S7. List of primers used in this study for quantitative PCR.
Genes included in 10 Gene Ontology (GO) categories closely related to drought stress. Red color indicates 48 genes differentially expressed only in BJ89 and BJ278 between drought stress and control treatments. Blue color indicates 49 genes differentially expressed between BJ89 and Nip or between BJ278 and Nip under drought stress treatments.
NAT pairs differentially expressed in Nip, BJ89, and BJ278 between drought stress and control treatments.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.