Identification and differential analysis of noncoding RNAs in response to drought in Phyllostachys aureosulcata f. spectabilis

The role of noncoding RNAs (ncRNAs) in plant resistance to abiotic stresses is increasingly being discovered. Drought stress is one of the most common stresses that affecting plant growth, and high intensity drought has a significant impact on the normal growth of plants. In this study, a high-throughput sequencing was performed on plant tissue samples of Phyllostachys aureosulcata f. spectabilis C. D. Chu et C. S. Chao by drought treatment for 0, 2, 4 and 6 days. The sequencing results were analysed bioinformatically. We detected 336,946 RNAs among all 12 samples, including 192,098 message RNAs (mRNAs), 142,761 long noncoding RNAs (lncRNAs), 1,670 circular RNAs (circRNAs), and 417 microRNAs (miRNAs). We detected 2,419 differentially expressed (DE) ncRNAs, including 213 DE circRNAs, 2,088 DE lncRNAs and 118 DE miRNAs. Then, we used Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) to functionally predict DE ncRNAs. The results showed that most DE ncRNAs are involved in the response to drought stress, mainly in biochemical reactions involved in some metabolites, as well as in organelle activities. In addition, we validated two random circRNAs and demonstrated their circularity. We also found a stable internal reference gene available for Phyllostachys aureosulcata f. spectabilis and validated the accuracy of this experiment by quantitative real-time polymerase chain reaction (qRT-PCR).

Introduction

Phyllostachys aureosulcata f. spectabilis C. D. Chu et C. S. Chao (P. aureosulcata) is a cultivated variant of Phyllostachys aureosulcata (Li et al., 2017). This kind of bamboo characterized by golden yellow poles, naturally long turquoise stripes. P. aureosulcata can stay green for all seasons, so it is especially suitable for cultivation under natural conditions in most areas (Zhang, 2014). According to the reasons above, the plant has become a rare bamboo treasure for landscaping and gardening worldwide. Especially after the Beijing Olympic Games, P. aureosulcata has been loved by people all over the world with its same name as the gold medal of the 2008 Beijing Summer Olympics and large area application in Olympic Park (Song et al., 2021b).

Drought is one of the environmental stresses that has the greatest worldwide impact on agricultural production (Kumar et al., 2018). In recent years, due to the drastic changes in global climate, the drought stresses encountered by plants under natural conditions have gradually increased in number and magnitude. Severe drought stress can cause physical damage and physiological or biochemical disorders, which in turn can alter plant morphology (Hussain et al., 2018), and this can be lethal to plants. Drought stress may also affect a range of physiological and biochemical reactions, including photosynthesis, chlorophyll synthesis, nutrient metabolism, ion uptake and transfer, respiratory action, and carbohydrate metabolism, inhibiting plant growth, and thus adversely affecting plant growth and yield (Li et al., 2011; Patel et al., 2020).

The plant response to drought stress is a complex biological process that involves a common dynamic change in metabolite composition and gene expression (Urano et al., 2009). Drought stress directly affects primary metabolism and leads to changes in the biosynthesis and transport of primary and secondary metabolites in plants (Ramakrishna and Ravishankar, 2011; Ma et al., 2020). At the gene expression level, plants have evolved various cascaded signalling networks that are used to regulate drought-responsive genes to produce various types of proteins, including transcription factors, enzymes, molecular chaperones, and other functional proteins, and ultimately achieve drought tolerance in plants (Hu and Xiong, 2014; Kumar et al., 2017; Kumar et al., 2021). Drought-responsive genes are involved in signalling cascades, transcriptional regulation, such as transcription factors and protein kinases/phosphatases, and the expression of functional proteins that protect cell membranes (Nakashima et al., 2014; Jogawat et al., 2021).

Noncoding RNAs (ncRNAs) are ubiquitous in plants, although for a long time, ncRNAs were viewed as pointless silent regions; later on, an increasing number of studies demonstrated that ncRNAs have a variety of regulatory effects (Bushati and Cohen, 2007; Bartel, 2009; Mercer et al., 2009). Identification and analysis in 15 diverse flowering plant species showed that 44% of ncRNAs have significant similarity to benchmark protein-coding or RNA genes and therefore have a high probability of being part of functional genes (Lloyd et al., 2018).

There are many types of ncRNAs, the main ones are microRNAs (miRNAs), long ncRNAs (lncRNAs), and circular RNAs (circRNAs), which usually interact with each other and commonly regulate the expression of target genes (Song et al., 2021a). ncRNAs have multiple functions. For example, ncRNAs are key regulators in several developmental processes such as leaf morphogenesis, vegetative phase change, flowering time and response to environmental cues, and the ncRNAs regulate the expression of abiotic stresses, such as drought, salt and high temperature, in plants facing different biotic stresses, salt, high temperature and other abiotic stresses, participation in plant immunity, and regulation in response to macronutrient stress in plants (Millar, 2020; Tahir Ul Qamar et al., 2020; Jin et al., 2021; Li et al., 2021).

In Bambusoideae, studies on ncRNAs have also made some progress. For example, they regulate the nutrients necessary for growth, influence the development of tissue and organs, and respond to abiotic stresses (Wang et al., 2016; Cheng et al., 2020; Wang et al., 2021), but no systematic studies about ncRNAs have been carried out in P. aureosulcata. Therefore, in this study, we systematically counted all types of ncRNAs in P. aureosulcata from the most common abiotic impact factor in growth: drought stress, focusing on ncRNAs involved in the stress response to drought stress, and performed KEGG/GO enrichment analysis. A stable reference gene was screened for P. aureosulcata and then verified by real-time quantitative PCR (qRT-PCR) using circRNAs as the point of penetration.

Materials and methods

Plant materials

Phyllostachys aureosulcata f. spectabilis C. D. Chu et C. S. Chao (P. aureosulcata) were grown in Dinghuo town, Yangzhou city, Jiangsu Province, China (119°38′10″E, 32°29′15″N). The average annual temperature is 14.8-15.3°C, and the light time is approximately 1896-2182 h per year. In addition, the annual precipitation can reach 1,048 mm, which is suitable for P. aureosulcata to grow. After digging the bamboo out of the soil, they were transplanted into planting pots (40 cm × 50 cm). After the seedlings were lowered, potted seedlings with good growth conditions were selected and treated with natural drought stress and divided into 4 groups with 3 replicates in each group. The drought treatment gradient was 0 days, 2 days, 4 days and 6 days. Leaves were collected from the same part of different plants, frozen with liquid nitrogen and saved at -80°C until further use. The soil water content is as follows: 0 day of drought treatment: 36.83%, 2 days of drought treatment: 29.22%, 4 days of drought treatment: 25.37%, 6 days of drought treatment: 22.15%.

RNA extraction, cDNA library construction and RNA sequencing

After dehydration with liquid nitrogen, the plant tissues were ground into powder, and total RNA was extracted using CLB, CTAB, and the RN40-EASYspin Plant RNA Rapid Extraction Kit from Aidlab Biotechnologies Co., Ltd., Beijing, China. Then, a Nanodrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA) was used for concentration testing, and an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) and LabChip GX Touch Nucleic Acid Analyser (PerkinElmer, Waltham, MA, USA) were used for integrity testing.

For cDNA library construction, a Ribo-off rRNA Depletion Kit (Plant) N409-02 was used. First, rRNA was removed from total RNA with rRNA probes and ferrite beads. Second, the interruption mixture was added to the library system reserved in step one to interrupt the rRNA-depleted RNA. Third, we synthesized the first strand of cDNA, synthesized the second chain according to the first chain and purified it. Fourth, end-point repair and dA-tailing of the above products were performed. Finally, reverse transcription was used to synthesize cDNA, followed by PCR amplification, polyacrylamide gel electrophoresis (PAGE) to separate the target DNA fragments, and cut gel recovery to obtain the cDNA library. While constructing the circRNA-seq libraries, add a step that removes the linear RNA with RNase R before the last step.

The quality of the constructed cDNA libraries was checked by the Qsep-400 method. The libraries that met the requirements were sequenced using Illumina NovaSeq 6000. The sequencing platform was an Illumina NovaSeq 6000 platform (San Diego, CA, USA), and the sequencing reagent was a NovaSeq 6000 S4 Reagent Kit (San Diego, CA, USA).

To identify the lncRNAs, three computational approaches including CPC2/CNCI/Pfam/CPAT were combined to sort non-protein coding RNA candidates from putative protein-coding RNAs in the unknown transcripts. Putative protein-coding RNAs were filtered out using a minimum length and exon number threshold. Transcripts with lengths more than 200 nt and have more than two exons been selected as lncRNA candidates (Li et al., 2015). To identify circRNAs, use CIRI software to compare with the reference gene sequences, generate SAM files, and scan PCC signals (paired chiastic clipping signals) from the SAM files for the CIGAR values analyzed in the SAM files. The CIGAR values in the junction read are characterized by xS/HyM or xMyS/H (Gao et al., 2015). Besides, to identify the sRNA. Raw reads of fastq format were firstly processed through in-house perl scripts. In this step, clean reads were obtained by removing reads containing adapter, reads containing ploy-N and low-quality reads from raw data. And reads were trimmed and cleaned by removing the sequences smaller than 18 nt or longer than 30 nt (Friedländer et al., 2012; Zhang et al., 2015).

Selection of internal reference genes and validation of circRNAs in P. aureosulcata

Moso bamboo (Phyllostachys edulis) is the most closely related and well-studied bamboo species of P. aureosulcata, so 11 candidate genes were selected for use as reference genes in moso bamboo (Fan et al., 2013; Qi et al., 2013; Wu et al., 2018). The conserved sequences were downloaded from the National Center for Biotechnology Information (NCBI) database (http://www.ncbi.nlm.nih.gov), and then Primer Premier 5 software was used to design primers for cloning (Supplementary Table 2). The annealing temperature of each primer was approximately 62°C, and the length was approximately 20-26 bp. Additionally, the length of each amplification product was between 100 bp and 180 bp. Respectively plucking the P. aureosulcata roots, stems, shoots, and leaves from different drought treatments in the Plant materials section (natural drought for 0 d, 2 d, 4 d, and 6 d), each tissue sample was obtained from at least three plants under good growth conditions. Different RNA was extracted by the TRIzol method (Total RNA Extraction Kit, TianMo Biotech, Tianjin, China), and 1 ng of each RNA was reverse transcribed to obtain cDNA for relative expression analysis of target genes by qRT-PCR. The total RNA was taken and directly reverse transcribed using PrimeScript RT Master Mix (Takara, Beijing, China) to obtain gDNA. Then, the total RNA was delinearized using Ribonuclease R (RNase R) and reverse transcribed using PrimeScript RT Master Mix to obtain cDNA with the linear RNA removed. A Nanodrop 2000 (Thermo Fisher Scientific) was used to check the RNA quality to ensure that the A260/A280 ratio of the RNA samples was between 1.9 and 2.1.

PCR amplification was carried out using cDNA as a template. The PCR products were visualized by agarose gel electrophoresis, and genes with clear and unique bands without primer-dimer formation were selected and confirmed by Sanger sequencing. After confirming that the sequences were correct, qRT-PCR was performed using 2×SYBR qPCR MasterMix (Real-Times Biotechnology, Beijing, China) on a Roche LightCycler 480 fluorescence quantification instrument (Basel, Switzerland). The qRT-PCR procedure was as follows: 95°C for 1 min, followed by 40 cycles of 95°C for 15 s and 60°C for 60 s. After obtaining the data, the results were imported into NormFinder software for statistical analysis, and the gene expression stability value (S value) was calculated (Andersen et al., 2004). The gene with the smallest S value was obtained as the most suitable internal reference gene for P. aureosulcata under the present experimental conditions.

The circRNAs were randomly selected, and the convergent and divergent primers were designed using Primer Premier 5. Then, cDNA and gDNA were used as templates with two sets of primers (Supplementary Table 2). The PCR products were also visualized by agarose gel electrophoresis and confirmed by Sanger sequencing. Using the system and apparatus above, qRT-PCR was performed with the selected PP2A as the reference gene, and the relative expression of circRNAs was calculated using the 2^-ΔΔCt method (Livak and Schmittgen, 2001).

Differential expression analysis

The DESeq R package (1.10.1) was used for differential expression analysis. DESeq provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P values were adjusted using Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with an adjusted P value <0.01 and absolute value of log₂(fold change) >1 found by DESeq were considered differentially expressed.

Functional prediction of DE RNAs

Gene Ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was implemented by the topGO R package. In addition, we used KEGG Orthology Based Annotation System (KOBAS) software to test the statistical enrichment of differentially expressed genes in KEGG pathways. Benjamini and Hochberg’s methods were used to adjust P values. In the detection of DE circRNAs, fold change greater than or equal to 2 and P values less than 0.05 were selected as screening criteria. The mRNA-miRNA-circRNA-lncRNA network was generated by Cytoscape (v3.7.2) (Shannon et al., 2003) and target gene prediction with TargetFinder software (Bo and Wang, 2005).

Results

Identification and characterization

To identify and associate ncRNAs with drought in P. aureosulcata, samples were taken from P. aureosulcata at 0 (P1), 2 (P2), 4 (P3) and 6 (P4) days of drought treatment, and RNA-Seq libraries were constructed, which were subjected to sequencing using the Illumina HiSeq 2500 platform. A total of 230.90 Gb of clean data was obtained from 12 samples, the percentage of Q30 bases in each sample was not less than 93.89%, and the guanine and cytosine (GC) content was approximately 44.58% (Supplementary Table 1).

Using the FASTQ data obtained from high-throughput sequencing, several RNAs were predicted and identified separately, and the quantities of several ncRNAs were counted in combination with known RNAs. The results are shown in Figure 1A. We detected 336,946 RNAs among all 12 samples, including 192,098 mRNAs, 142,761 lncRNAs, 1,670 circRNAs, and 417 miRNAs. Then, we used 3 different software programs to analyse different kinds of ncRNAs to obtain more details. CircRNAs were predicted using circRNA identifier (CIRI) (Gao et al., 2015), and the distribution statistics of the length of predicted circRNAs were carried out. As shown in Figure 1B, most of the circRNAs were between 200-800 bp in length. The type statistics of circRNAs performed in Figure 1C also showed that most of the circRNAs were classified as exon circRNAs. Then, using HISAT2 (Li et al., 2015) to predict lncRNAs. The predicted lncRNAs were distributed with length and class statistics, as shown in Figures 1D, E. Most of the lncRNAs were between 400-800 bp in length and were classified as exon lncRNAs. Finally, using miRDeep2 software, adjusting and changing its parameters and scoring system to make it suitable for the prediction of plant miRNAs (Friedländer et al., 2012; Zhang et al., 2015). As shown in Figure 1F, most of the miRNAs were 21 or 24 bp in length.

FIGURE 1

Figure 1 Characterization of RNAs. (A) The number of each type of RNA in all samples. (B) The length distribution of circRNAs. (C) The type and percentage of circRNAs in each sample. (D) The type and percentage of lncRNAs in each sample. (E) The length distribution of lncRNAs. (F) The length distribution of miRNAs in each sample. Drought treatment with three biological repetitions each: 0 day (P11, P12, P13), 2 days of drought treatment (P21, P22, P23), 4 days of drought treatment (P31, P32, P33) and 6 days of drought treatment (P41, P42, P43).

Differentially expressed ncRNAs under drought treatment

A total of 2,419 differentially expressed ncRNAs were detected under different gradients of drought treatment. Among them, a total of 213 DE circRNAs were detected. There were also 2,088 DE lncRNAs. The total number of DE miRNAs with differential epistatic responses was 118 (Figure 2). The clustering analysis of differentially expressed lncRNAs, circRNAs and miRNAs was performed, and the results are shown in Figure 2.

FIGURE 2

Figure 2 Differential expression analysis of ncRNAs in response to drought stress treatment. Drought treatment with three biological repetitions each: 0 day (P11, P12, P13), 2 days of drought treatment (P21, P22, P23), 4 days of drought treatment (P31, P32, P33) and 6 days of drought treatment (P41, P42, P43). (A) The regulation of different kinds of DERNA between each treatment. (B) Venn diagram of DE circRNAs at six group of treatment. (C) Heatmap of expression of all DE circRNAs in all samples, the horizontal coordinates in the diagram represent the sample names and the clustering results of the samples, the vertical coordinates represent the differential genes and the clustering results of the genes. Different columns in the graph represent different samples, and different rows represent different genes. The color represents the expression level of the gene in the sample log10(FPKM+0.000001). (D) Venn diagram of DE circRNAs at six group of treatment. (E) Heatmap of expression of all DE lncRNAs in all samples. (F) Venn diagram of DE miRNAs at six group of treatment. (G) Heatmap of expression of all DE miRNAs in all samples.

Reference gene selection and quantitative real-time polymerase chain reaction

Among the 11 candidate reference genes, 6 candidate genes were measured by qRT-PCR products of agarose gel electrophoresis and showed a single specific band with the same size of the target fragment and no primer dimer below the band. The six selected reference genes were Phosphoprotein Phosphatase 2A (PP2A), Actin2-1 (ACT2-1), Elongation Factor 1 Alpha (EF1α), Ubiquitin (UBQ), Actin-1 (ACT1), and Clathrin Adaptor Complexes Medium Subunit (CAC) (Figure 3A). Their primers were amplified with high efficiency and met the criteria of qRT-PCR experiments for gene expression stability analysis. For the six selected reference genes, qRT-PCR analysis was performed with cDNAs under different drought treatments (0 days, 2 days, 4 days and 6 days) and cDNAs from different plant organs (roots, stems, leaves and shoots) as templates. The results were imported into NormFinder software for statistical analysis (Andersen et al., 2004), and the gene expression stability values (S values) were calculated (Table 1). The gene with the smallest S value under both situations was obtained as PP2A, and PP2A was identified as the most suitable internal reference gene for P. aureosulcata under the present experimental conditions.

FIGURE 3

Figure 3 Reference gene selection and qRT-PCR validation of four DE circRNAs (A) PCR amplification of 11 reference genes. (B-E) Relative expression of four selected DE circRNAs measured by RT-qPCR. For RT-qPCR, three replicates were performed. Asterisks on bars indicate significant differences between different stages using a t-test. Error bars represent the standard error of the mean. Drought treatment: 0 day (P1), 2 days of drought treatment (P2), 4 days of drought treatment (P3) and 6 days of drought treatment (P4). * p < 0.05, indicates a significant difference. ** p < 0.01, indicates a profoundly significant difference. *** p < 0.001, indicates a profoundly significant difference.

TABLE 1

Table 1 Ranking order of reference genes determined by NormFinder under drought treatments.

In addition, to verify the reliability of the gene expression profile, we randomly selected four DE circRNAs for expression analysis (Figure 3). qRT-PCR results showed expression patterns consistent with the RNA-seq results.

Validation of circRNAs

To verify the accuracy of our identified circRNAs, we randomly selected two circRNAs using PCR and Sanger sequencing to validate them and successfully confirmed their circularity. Unlike linear RNAs, circRNAs are a loop, so two sets of primers (convergent and divergent) and two templates (genomic DNA [gDNA] and complementary DNA [cDNA]) were used to verify it (Figure 4A). The normal convergent primers can obtain amplification bands in both gDNA and cDNA, while the reverse divergent primers can only obtain amplification bands in cDNA (Figure 4).

FIGURE 4

Figure 4 Validation of circRNAs by PCR and Sanger sequencing. (A) Schematic representation of the validation of circRNA1: hic_scaffold_22:10652336|10653429. circRNA1 is formed after the circularization of 4 exons in the PH02Gene45872 gene. (B) Schematic representation of the validation of circRNA2: hic_scaffold_8:14927508|14929885. circRNA2 is formed after the circularization of 2 exons in the PH02Gene31681 gene. PCR amplification was performed after designing two sets of primers (convergent and divergent) for these two genes; clear bands could be observed in the products by agarose gel electrophoresis (AGE), and the amplified products could be sequenced to obtain the complete sequences at the first and last junctions.

Functional prediction of DE ncRNAs by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes

To further understand the functions and characteristics of DE ncRNAs, enrichment analyses of all DE ncRNAs, DE circRNAs, DE lncRNAs, and DE miRNAs were performed using the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (Figures 5, 6 and Supplementary Tables 3, 4).

FIGURE 5

Figure 5 GO and KEGG enrichment analyses of the host genes of all DE RNAs. (A) GO enrichment analysis of the host genes of all DE ncRNAs. The horizontal coordinate is the GO classification, the left side of the vertical coordinate is the percentage of the number of genes enriched in that GO term, and the right side is the number of genes. (B) KEGG enrichment analysis of the host genes of all DE ncRNAs. Each circle in Figure 5 represents a KEGG pathway, the vertical coordinate indicates the top 20 pathway names, and the horizontal coordinate is the enrichment factor, which represents the ratio of the proportion of differential genes annotated to a pathway (gene ratio) to the proportion of all genes annotated to that pathway (background ratio). The colour of the circle represents the q value, which is the P value after correction for multiple hypothesis testing, and the size of the circle indicates the number of genes enriched in the pathway.

FIGURE 6

Figure 6 GO and KEGG enrichment analyses of the host genes of DE RNAs separately. (A) GO enrichment analysis of the host genes of DE circRNAs. The horizontal coordinate is the GO classification, the left side of the vertical coordinate is the percentage of the number of genes enriched in that GO term, and the right side is the number of genes. (B) KEGG enrichment analysis of the host genes of DE circRNAs. Each circle in Figure 6 represents a KEGG pathway, the vertical coordinate indicates the top 20 pathway names, and the horizontal coordinate is the enrichment factor, which represents the ratio of the proportion of differential genes annotated to a pathway (gene ratio) to the proportion of all genes annotated to that pathway (background ratio). The colour of the circle represents the q value, which is the P value after correction for multiple hypothesis testing, and the size of the circle indicates the number of genes enriched in the pathway. (C) GO enrichment analysis of the host genes of DE lncRNAs. (D) KEGG enrichment analysis of the host genes of DE circRNAs. (E) GO enrichment analysis of the host genes of DE miRNAs. (F) KEGG enrichment analysis of the host genes of DE miRNAs.

According to the results of GO enrichment analysis, the host genes were involved in three categories: biological process (BP), cellular component (CC), and molecular function (MF). In the BP section, most host genes were enriched in metabolic process (GO: 0008152), cellular process (GO: 0044763), single organism process (GO: 0006465) and biological regulation (GO: 0065007). In the CC category, most host genes were enriched in cell (GO: 0005623), cell part (GO: 0044464), organelle (GO: 0043226) and membrane (GO: 0016020). In addition, in the MF category, most host genes were enriched in binding (GO: 0005488), catalytic activity (GO: 0003824), transporter activity (GO: 0005215) and nucleic acid binding transcription factor activity (GO: 0003700).

For the KEGG enrichment analysis, the most closely related metabolic pathways of the host genes of all DE ncRNAs were glycolysis/gluconeogenesis (ko00010), endocytosis (ko04144), Carotenoid biosynthesis (ko00906), Terpenoid backbone biosynthesis (ko00900), Phosphonate and phosphinate metabolism (ko00440). Besides, the results showed that the most closely related metabolic pathways of the host genes of DE circRNAs under drought treatment were fatty acid biosynthesis (ko00061), fatty acid metabolism (ko01212), peroxisome (ko04146), fatty acid degradation (ko00071) and RNA transport (ko03013). In addition, the most closely related metabolic pathways of the host gene of DE lncRNAs are phenylpropanoid biosynthesis (ko00940), DNA replication (ko03030), starch and sucrose metabolism (ko00500), pentose and glucuronate interconversions (ko00040) and plant hormone signal transduction (ko04075). Meanwhile, the most closely related metabolic pathways of the host genes of DE miRNAs are steroid biosynthesis (ko00905), biosynthesis of various secondary metabolites-Part 2 (ko00998), inositol phosphate metabolism (ko00562), flavonoid biosynthesis (ko00941) and photosynthesis-antenna proteins (ko00196).

Functional prediction of DE ncRNAs based on the mRNA-miRNA-circRNA-lncRNA network

To further explore the functions of DE ncRNAs, we first predicted DE circRNAs and DE lncRNAs as miRNA decoys and mRNAs as miRNA targets. Then, based on the potential relationship between these RNAs, three mRNA-miRNA-circRNA-lncRNA regulatory networks were constructed by screening the differentially expressed mRNAs (Figure 7). The reciprocal relationship networks of P1 treatment with drought treatment P2 (Supplementary Figure 1), P3 (Supplementary Figure 1) and P4 (Figure 7). Meanwhile, the mRNAs involved in the above three mRNA-miRNA-circRNA-lncRNA regulatory networks were annotated and analysed. We found that these mRNAs were mainly annotated to the following GO entries:

FIGURE 7

Figure 7 The mRNA-miRNA-circRNA-lncRNA interaction network of P1 treatment with drought treatment P4. Yellow nodes: miRNAs. Red nodes: circRNAs that may be miRNA decoys. Blue nodes: lncRNAs that may be miRNA decoys. Pink nodes: mRNAs that may be miRNA targets.

Molecular Function: DNA binding (GO: 0003677), nucleic acid binding (GO: 0003676), zinc ion binding (GO: 0008270), RNA binding (GO: 0003723), RNA-directed DNA polymerase activity (GO: 0003964), aspartic-type endopeptidase activity (GO: 0004190), transporter activity (GO: 0005215), phenylalanine ammonia-lyase activity (GO: 0045548). Cellular Component: Nucleus (GO: 0005634), Plastid (GO: 0009536), Integral component of membrane (GO: 0016021), Membrane (GO: 0016020), Biological Process: regulation of transcription, DNA-templated (GO: 0006355), Mitochondrion (GO: 0005739), Cytoplasm (GO: 0005737). Biological Process: DNA integration (GO: 0015074), Regulation of transcription, DNA-templated (GO: 0006355), Auxin-activated signalling pathway (GO: 0009734), RNA-dependent DNA biosynthetic process (GO: 0006278), DNA recombination (GO: 0006310), proteolysis (GO: 0006508), L-phenylalanine catabolic process (GO: 0006559), Cinnamic acid biosynthetic process (GO: 0009800).

Based on the above two chapters, the KEGG/GO enrichment analysis showed that non-coding RNAs are mainly involved in the following biological processes: 1, metabolism of primary/secondary metabolites. This is the class with the highest number of enriched ncRNAs and involving GO/KEGG entries. By regulating the synthesis and metabolism of some metabolites closely related to the daily life activities of plants and key compounds in key biological cycles (e.g., saccharides, alcohols, proteins, lipids, etc.) to counteract the effects of external environmental water changes. 2, Organelle activity. By regulating the activity of organelles such as peroxisome, plastid, and mitochondrion, the strength of relevant biochemical reactions in plants is mobilized to increase/decrease, thus enhancing the resistance of plants to external water changes. 3, hormones. By regulating plant hormone signal transduction to make changes in plant physiological and biochemical metabolism, thus improving drought resistance and reducing plant damage. 4, Increase drought resistance by regulating RNA and DNA synthesis/recombination/translation, etc. 5, Regulate the ions in the tissues. By regulating the activity of cell membrane, as well as regulating the state of ions such as Zn²⁺ and Ca²⁺ in the cell, and regulating the osmotic pressure in the tissues to resist external water changes.

Discussion

As global warming increases, environmental temperatures are gradually rising, which greatly increases the probability and duration of drought. As a common greenery bamboo, P. aureosulcata is planted mostly in urban environments and grows naturally, which further deepens the impact of drought stress on it. Meanwhile, the role of ncRNAs has been gaining attention in recent decades, from inconsequential transcriptional “noise” to “a new continent in the RNA world”. ncRNAs are mainly classified into long noncoding RNAs (lncRNAs), which are longer than 200 nt, and small noncoding RNAs (sRNAs), which are smaller than 200 nt (Gelaw and Sanan-Mishra, 2021). Many experimental reports have shown that ncRNAs and their target genes play important roles in plant development, environmental responses, and biotic and abiotic stress responses (Ohtani, 2017).

With the advancement of high-throughput sequencing technology, many new ncRNA transcripts have been identified in different species (Di et al., 2014; Wang et al., 2015; Deng et al., 2018; Yuan et al., 2018). Bioinformatics analysis has proven to be an important tool for exploring differences in plant responses to drought stress and has been applied to several species, including Arabidopsis (Qin et al., 2017), rice (Mutum et al., 2016), maize (Zhang et al., 2014), tomato (Zhong et al., 2013), sorghum (Fracasso et al., 2016), and coffee (Mofatto et al., 2016), and many drought-related DE ncRNAs have been identified. These studies revealed the complexity of the regulation between drought stress and ncRNAs and revealed that ncRNAs play important roles in many important biological processes. Therefore, understanding the regulatory mechanisms of ncRNAs in response to drought will provide a molecular basis for plant resistance studies. However, genomic identification and characterization of known and novel ncRNAs under drought stress of P. aureosulcata are still lacking. Therefore, in this research, we conducted an overall statistical analysis, identification and analysis of all species of DE ncRNAs in P. aureosulcata under drought treatment and focused on analysing circRNAs in addition to counting conventional lncRNAs and miRNAs. A total of 2,419 differentially expressed ncRNAs were detected under different gradients of drought treatment, including 213 DE circRNAs, 2,088 DE lncRNAs and 118 DE miRNAs. Among these DE ncRNAs, the number of upregulated versus downregulated DE ncRNAs was largely maintained at 1:1, which is consistent with previous studies. The DE ncRNAs identified in this study originated from exons, introns, and intergenic regions, with most DE ncRNAs, especially DE circRNAs, originating from a single exon, which may be related to the current mechanism of circRNA formation in plants: exon skipping events (Conn et al., 2017). In addition, the expression of many circRNAs was extremely low, a feature that may be an essential feature common to circRNAs in plants (Zhang et al., 2020).

In many plant species, ncRNAs have been found to be involved in regulating the expression of their host genes, regulating the physiological and biochemical activities of the plant, and thus regulating the growth life of the plant. In the face of drought, ncRNAs perform different functions to regulate the state of the plant itself. For example, in Arabidopsis thaliana, drought-induced lncRNA (DRIR) was confirmed as a novel positive regulator of the plant response to drought and salt stresses (Qin et al., 2017), and miR165/166 use abscisic acid (ABA) signalling to modify plant xylem morphology under conditions of environmental stress (Ramachandran et al., 2018). In cassava, lncRNAs were confirmed to play crucial roles in MT-mediated drought stress responses (Ding et al., 2019). In maize, DRIR represses ZmNAC111 expression and enhances drought tolerance through RNA-directed DNA methylation (Pang et al., 2019), and several microRNAs can regulate phosphate uptake and affect the growth of primary roots in response to nutrient deficiencies in response to water-deficit stress (Seeve et al., 2019). Therefore, in this study, GO enrichment analysis and KEGG enrichment analysis were performed on host genes of DERNAs to investigate the function of DERNAs in the drought response. GO analysis showed that the host genes were enriched mainly in metabolic process (GO: 0008152), cellular process (GO: 0044763), cell (GO: 0005623), binding (GO: 0005488), and catalytic activity (GO: 0003824). KEGG analysis showed that the most closely related metabolic pathways of the host genes of DERNAs under drought treatment were fatty acid biosynthesis (ko00061), fatty acid metabolism (ko01212), peroxisome (ko04146), fatty acid degradation (ko00071), phenylpropanoid biosynthesis (ko00940), DNA replication (ko03030), starch and sucrose metabolism (ko00500), pentose and glucuronate interconversions (ko00040), plant hormone signal transduction (ko04075), steroid biosynthesis (ko00905) and biosynthesis of various secondary metabolites-Part 2 (ko00998). These pathways are clearly related to the synthesis and degradation of primary and secondary metabolites and involve synergistic interactions between multiple organelles and hormonal signalling. Therefore, we hypothesized that these ncRNAs regulate the functions of various organelles and regulate the levels of various hormones, leading to changes in the metabolism of multiple metabolites and thus resisting the lack of water within the environment.

Data availability statement

The original contributions presented in the study are publicly available. This data can be found here: NCBI, PRJNA885522 and PRJNA885524.

Author contributions

XL and YY planned this experiment. YY, YG and YL performed the experiments. YY and YG analyzed the data. YY write this article. All authors contributed to the article and approved the submitted version.

Funding

This research was funded by the National Key R&D Program of China, 2021YFD2200504_4.

Conflict of interest

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

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.1040470/full#supplementary-material

References

Andersen, C. L., Jensen, J. L., Ørntoft, T. F. (2004). Normalization of real-time quantitative reverse transcription-PCR data: A model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 64 (15), 5245–5250. doi: 10.1158/0008-5472.CAN-04-0496

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartel, D. P. (2009). MicroRNAs: target recognition and regulatory functions. Cell 136 (2), 215–233. doi: 10.1016/j.cell.2009.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Bo, X., Wang, S. (2005). TargetFinder: A software for antisense oligonucleotide target site selection based on MAST and secondary structures of target mRNA. Bioinformatics 21 (8), 1401–1402. doi: 10.1093/bioinformatics/bti211

PubMed Abstract | CrossRef Full Text | Google Scholar

Bushati, N., Cohen, S. M. (2007). microRNA functions. Annu. Rev. Cell Dev. Biol. 23, 175–205. doi: 10.1146/annurev.cellbio.23.090506.123406

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, Z., Hou, D., Ge, W., Li, X., Xie, L., Zheng, H., et al. (2020). Integrated mRNA, MicroRNA transcriptome and degradome analyses provide insights into stamen development in moso bamboo. Plant Cell Physiol. 61 (1), 76–87. doi: 10.1093/pcp/pcz179

PubMed Abstract | CrossRef Full Text | Google Scholar

Conn, V. M., Hugouvieux, V., Nayak, A., Conos, S. A., Capovilla, G., Cildir, G., et al. (2017). A circRNA from SEPALLATA3 regulates splicing of its cognate mRNA through r-loop formation. Nat. Plants 3(17053), 12017. doi: 10.1038/nplants.2017.53

CrossRef Full Text | Google Scholar

Deng, F., Zhang, X., Wang, W., Yuan, R., Shen, F. (2018). Identification of Gossypium hirsutum long non-coding RNAs (lncRNAs) under salt stress. BMC Plant Biol. 18 (1), 23. doi: 10.1186/s12870-018-1238-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, Z., Wu, C., Tie, W., Yan, Y., He, G., Hu, W. (2019). Strand-specific RNA-seq based identification and functional prediction of lncRNAs in response to melatonin and simulated drought stresses in cassava. Plant Physiol. Biochem. 140, 96–104. doi: 10.1016/j.plaphy.2019.05.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Di, C., Yuan, J., Wu, Y., Li, J., Lin, H., Hu, L., et al. (2014). Characterization of stress-responsive lncRNAs in Arabidopsis thaliana by integrating expression, epigenetic and structural features. Plant J. 80 (5), 848–861. doi: 10.1111/tpj.12679

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, C., Ma, J., Guo, Q., Li, X., Wang, H., Lu, M. (2013). Selection of reference genes for quantitative real-time PCR in bamboo (Phyllostachys edulis). PloS One 8 (2), e56573. doi: 10.1371/journal.pone.0056573

PubMed Abstract | CrossRef Full Text | Google Scholar

Fracasso, A., Trindade, L. M., Amaducci, S. (2016). Drought stress tolerance strategies revealed by RNA-seq in two sorghum genotypes with contrasting WUE. BMC Plant Biol. 16 (1), 115. doi: 10.1186/s12870-016-0800-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedländer, M. R., Mackowiak, S. D., Li, N., Chen, W., Rajewsky, N. (2012). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40 (1), 37–52. doi: 10.1093/nar/gkr688

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Y., Wang, J., Zhao, F. (2015). CIRI: An efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 16, 4. doi: 10.1186/s13059-014-0571-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Gelaw, T. A., Sanan-Mishra, N. (2021). Non-coding RNAs in response to drought stress. Int. J. Mol. Sci. 22 (22), 12519. doi: 10.3390/ijms222212519

PubMed Abstract | CrossRef Full Text | Google Scholar

Hussain, H. A., Hussain, S., Khaliq, A., Ashraf, U., Anjum, S. A., Men, S., et al. (2018). Chilling and drought stresses in crop plants: Implications, cross talk, and potential management opportunities. Front. Plant Sci. 9. doi: 10.3389/fpls.2018.00393

CrossRef Full Text | Google Scholar

Hu, H., Xiong, L. (2014). Genetic engineering and breeding of drought-resistant crops. Annu. Rev. Plant Biol. 65, 715–741. doi: 10.1146/annurev-arplant-050213-040000

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, J., Lu, P., Xu, Y., Li, Z., Yu, S., Liu, J., et al. (2021). PLncDB V2.0: A comprehensive encyclopedia of plant long noncoding RNAs. Nucleic Acids Res. 49 (D1), D1489–D1495. doi: 10.1093/nar/gkaa910

PubMed Abstract | CrossRef Full Text | Google Scholar

Jogawat, A., Yadav, B., Chhaya., Lakra, N., Singh, A. K., Narayan, O. P. (2021). Crosstalk between phytohormones and secondary metabolites in the drought stress tolerance of crop plants: A review. Physiol. Plant 172 (2), 1106–1132. doi: 10.1111/ppl.13328

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, R., Bohra, A., Pandey, A. K., Pandey, M. K., Kumar, A. (2017). Metabolomics for plant improvement: Status and prospects. Front. Plant Sci. 8. doi: 10.3389/fpls.2017.01302

CrossRef Full Text | Google Scholar

Kumar, M., Kumar Patel, M., Kumar, N., Bajpai, A. B., Siddique, K. H. M. (2021). Metabolomics and molecular approaches reveal drought stress tolerance in plants. Int. J. Mol. Sci. 22 (17), 9108. doi: 10.3390/ijms22179108

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, M., Yusuf, M. A., Nigam, M., Kumar, M. (2018). An update on genetic modification of chickpea for increased yield and stress tolerance. Mol. Biotechnol. 60 (8), 651–663. doi: 10.1007/s12033-018-0096-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Gao, J., Fan, R. (2017). Changes in physiological indices and leaf structure of Phyllostachys aureosulacata f. spectabilis, Ph.vivax f. aureocaulis, Ph. vivax f. huangwenzhu during winter in Beijing. J. Fujian Agric. Sci. Tech. 46 (5), 527–533. doi: 10.13323/j.cnki.j.fafu(nat.sci.).2017.05.008

CrossRef Full Text | Google Scholar

Li, C., Jiang, D., Wollenweber, B., Li, Y., Dai, T., Cao, W. (2011). Waterlogging pretreatment during vegetative growth improves tolerance to waterlogging after anthesis in wheat. Plant Sci. 180 (5), 672–678. doi: 10.1016/j.plantsci.2011.01.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Ma, W., Zeng, P., Wang, J., Geng, B., Yang, J., et al. (2015). LncTar: a tool for predicting the RNA targets of long noncoding RNAs. Brief. Bioinform. 16 (5), 806–812. doi: 10.1093/bib/bbu048

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Tian, P., Huang, T., Huang, J. (2021). Noncoding-RNA-Mediated regulation in response to macronutrient stress in plants. Int. J. Mol. Sci. 22 (20), 11205. doi: 10.3390/ijms222011205

PubMed Abstract | CrossRef Full Text | Google Scholar

Livak, K. J., Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2(-delta delta C(T)) method. Methods 25 (4), 402–408. doi: 10.1006/meth.2001.1262

PubMed Abstract | CrossRef Full Text | Google Scholar

Lloyd, J. P., Tsai, Z. T., Sowers, R. P., Panchy, N. L., Shiu, S. H. (2018). A model-based approach for identifying functional intergenic transcribed regions and noncoding RNAs. Mol. Biol. Evol. 35 (6), 1422–1436. doi: 10.1093/molbev/msy035

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, X. J., Yu, T. F., Li, X. H., Cao, X. Y., Ma, J., Chen, J., et al. (2020). Overexpression of GmNFYA5 confers drought tolerance to transgenic arabidopsis and soybean plants. BMC Plant Biol. 20 (1), 123. doi: 10.1186/s12870-020-02337-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Mercer, T. R., Dinger, M. E., Mattick, J. S. (2009). Long non-coding RNAs: insights into functions. Nat. Rev. Genet. 10 (3), 155–159. doi: 10.1038/nrg2521

PubMed Abstract | CrossRef Full Text | Google Scholar

Millar, A. A. (2020). The function of miRNAs in plants. Plants (Basel) 9 (2), 198. doi: 10.3390/plants9020198

CrossRef Full Text | Google Scholar

Mofatto, L. S., Carneiro Fde, A., Vieira, N. G., Duarte, K. E., Vidal, R. O., Alekcevetch, J. C., et al. (2016). Identification of candidate genes for drought tolerance in coffee by high-throughput sequencing in the shoot apex of different Coffea arabica cultivars. BMC Plant Biol. 16, 94. doi: 10.1186/s12870-016-0777-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Mutum, R. D., Kumar, S., Balyan, S., Kansal, S., Mathur, S., Raghuvanshi, S. (2016). Identification of novel miRNAs from drought tolerant rice variety nagina 22. Sci. Rep. 6, 30786. doi: 10.1038/srep30786

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakashima, K., Yamaguchi-Shinozaki, K., Shinozaki, K. (2014). The transcriptional regulatory network in the drought response and its crosstalk in abiotic stress responses including drought, cold, and heat. Front. Plant Sci. 5. doi: 10.3389/fpls.2014.00170

PubMed Abstract | CrossRef Full Text | Google Scholar

Ohtani, M. (2017). Expanding the plant non-coding RNA world. J. Plant Res. 130 (1), 3–5. doi: 10.1007/s10265-016-0896-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Pang, J., Zhang, X., Ma, X., Zhao, J. (2019). Spatio-temporal transcriptional dynamics of maize long non-coding RNAs responsive to drought stress. Genes (Basel) 10 (2), 138. doi: 10.3390/genes10020138

CrossRef Full Text | Google Scholar

Patel, M. K., Kumar, M., Li, W., Luo, Y., Burritt, D. J., Alkan, N., et al. (2020). Enhancing salt tolerance of plants: From metabolic reprogramming to exogenous chemical treatments and molecular approaches. Cells 9 (11), 2492. doi: 10.3390/cells9112492

CrossRef Full Text | Google Scholar

Qi, F., Hu, T., Peng, Z., Gao, J. (2013). Screening of reference genes used in qRT-PCR and expression analysis of PheTFL1 gese in moso bamboo. Acta Bot. Boreal. -Occident. Sin. 33 (1), 48–52. doi: 10.3969/j.issn.1000-4025.2013.01.010

CrossRef Full Text | Google Scholar

Qin, T., Zhao, H., Cui, P., Albesher, N., Xiong, L. (2017). A nucleus-localized long non-coding RNA enhances drought and salt stress tolerance. Plant Physiol. 175 (3), 1321–1336. doi: 10.1104/pp.17.00574

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramachandran, P., Wang, G., Augstein, F., de Vries, J., Carlsbecker, A. (2018). Continuous root xylem formation and vascular acclimation to water deficit involves endodermal ABA signalling via miR165. Development 145 (3), dev.159202. doi: 10.1242/dev.159202

CrossRef Full Text | Google Scholar

Ramakrishna, A., Ravishankar, G. A. (2011). Influence of abiotic stress signals on secondary metabolites in plants. Plant Signal. Behav. 6 (11), 1720–1731. doi: 10.4161/psb.6.11.17613

PubMed Abstract | CrossRef Full Text | Google Scholar

Seeve, C. M., Sunkar, R., Zheng, Y., Liu, L., Liu, Z., McMullen, M., et al. (2019). Water-deficit responsive microRNAs in the primary root growth zone of maize. BMC Plant Biol. 19 (1), 447. doi: 10.1186/s12870-019-2037-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13 (11), 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, L., Fang, Y., Chen, L., Wang, J., Chen, X. (2021a). Role of non-coding RNAs in plant immunity. Plant Commun. 2 (3), 100180. doi: 10.1016/j.xplc.2021.100180

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, H., Zhang, A., Wang, R., Qi, L. (2021b). Introduction history, cultivation research and landscape application of bamboo in Beijing. World Forestry Res. 34 (6), 67–71. doi: 10.13348/j.cnki.sjlyyj.2021.0061.y

CrossRef Full Text | Google Scholar

Tahir Ul Qamar, M., Zhu, X., Khan, M. S., Xing, F., Chen, L. L. (2020). Pan-genome: A promising resource for noncoding RNA discovery in plants. Plant Genome 13 (3), e20046. doi: 10.1002/tpg2.20046

PubMed Abstract | CrossRef Full Text | Google Scholar

Urano, K., Maruyama, K., Ogata, Y., Morishita, Y., Takeda, M., Sakurai, N., et al. (2009). Characterization of the ABA-regulated global responses to dehydration in arabidopsis by metabolomics. Plant J. 57 (6), 1065–1078. doi: 10.1111/j.1365-313X.2008.03748.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Hou, Y., Wang, Y., Zhao, H. (2021). Integrative lncRNA landscape reveals lncRNA-coding gene networks in the secondary cell wall biosynthesis pathway of moso bamboo (Phyllostachys edulis). BMC Genomics 22 (1), 638. doi: 10.1186/s12864-021-07953-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, T. Z., Liu, M., Zhao, M. G., Chen, R., Zhang, W. H. (2015). Identification and characterization of long non-coding RNAs involved in osmotic and salt stress in Medicago truncatula using genome-wide high-throughput sequencing. BMC Plant Biol. 15, 131. doi: 10.1186/s12870-015-0530-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Zhao, H., Chen, D., Li, L., Sun, H., Lou, Y., et al. (2016). Characterization and primary functional analysis of a bamboo NAC gene targeted by miR164b. Plant Cell Rep. 35 (6), 1371–1383. doi: 10.1007/s00299-016-1970-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, L., Zhang, Z., Zhu, F., Zhang, M., Yu, Y., Qiu, Y., et al. (2018). The selection of endogenous reference genes in Phyllostachys edulis for qRT-PCR analysis. J. Agric. Biotechnol. 26 (3), 502–510. doi: 10.3969/j.issn.1674-7968.2018.03.016

CrossRef Full Text | Google Scholar

Yuan, J., Li, J., Yang, Y., Tan, C., Zhu, Y., Hu, L., et al. (2018). Stress-responsive regulation of long non-coding RNA polyadenylation in Oryza sativa. Plant J. 93 (5), 814–827. doi: 10.1111/tpj.13804

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y. (2014). Study on shooting and bamboo of Phyllostachys aureosulcata f. spectabilis introduced into south area of henan province. Jour. Fujian Forestry Sci. Tech. 41 (4), 26–29. doi: 10.13428/j.cnki.fjlk.2014.04.006

CrossRef Full Text | Google Scholar

Zhang, W., Han, Z., Guo, Q., Liu, Y., Zheng, Y., Wu, F., et al. (2014). Identification of maize long non-coding RNAs responsive to drought stress. PloS One 9 (6), e98958. doi: 10.1371/journal.pone.0098958

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z., Jiang, L., Wang, J., Gu, P., Chen, M. (2015). MTide: an integrated tool for the identification of miRNA-target interaction in plants. Bioinformatics 31 (2), 290–291. doi: 10.1093/bioinformatics/btu633

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Liu, R., Zhu, Y., Gong, J., Yin, S., Sun, P., et al. (2020). Identification and characterization of circRNAs responsive to methyl jasmonate in Arabidopsis thaliana. Int. J. Mol. Sci. 21 (3), 792. doi: 10.3390/ijms21030792

CrossRef Full Text | Google Scholar

Zhong, S., Fei, Z., Chen, Y. R., Zheng, Y., Huang, M., Vrebalov, J., et al. (2013). Single-base resolution methylomes of tomato fruit development reveal epigenome modifications associated with ripening. Nat. Biotechnol. 31 (2), 154–159. doi: 10.1038/nbt.2462

PubMed Abstract | CrossRef Full Text | Google Scholar