Integration of mRNA and miRNA analysis reveals the differentially regulatory network in two different Camellia oleifera cultivars under drought stress

Camellia oleifera Abel. (C. oleifera) is an edible oil tree species that provide an important guarantee for targeted poverty alleviation strategy in China. Severe difficulties in irrigation leading to drought stress have become a major obstacle to the development of the C. oleifera planting industry. Breeding of drought-tolerant cultivars is the main idea for solving the problem of water shortage stress in C. oleifera cultivation. The photosynthetic physiology traits of C. oleifera cultivars ‘Xianglin No.1’ and ‘Hengdong No.2’ were affected by drought stress to different degrees, which demonstrated that the two cultivars suffered different degrees of damage. In the present study, we applied mRNA-seq and miRNA-seq to analyze the difference in molecular responses between drought stress and control, drought-tolerant and -sensitive cultivars, at mRNA and miRNA levels. The differentially expressed genes (DEGs) involved in photosynthesis-related, porphyrin, and chlorophyll metabolism, circadian rhythm system, and plant hormone signal transduction pathways were identified that might be candidates for drought stress tolerance genes. Subsequently, the miRNA-mRNA regulatory networks connected the differentially expressed miRNAs (DEMs) to their predicted target genes were established. miR398 and miR408-3p in C. oleifera showed that associated with the response to drought stress by negatively regulating genes encoding Downy Mildew Resistance 6 (DMR6) and Enhanced Disease Resistance 2 (EDR2), respectively, which might further improve drought tolerance via crosstalk between different stress-responsive pathways. The quantitation results of miRNA and mRNA were validated by quantitative real-time PCR (qRT-PCR). In summary, the integrated mRNA-seq and miRNA-seq analysis deepen our understanding of the regulatory network response to drought stress and variety-specific responses improving drought tolerance in C. oleifera.

Introduction

Camellia oleifera Abel. (C. oleifera) is an edible oil tree species with growing commercial, medic, cosmetic and ornamental values in recent years that provide an important guarantee for targeted poverty alleviation strategy in China (He et al., 2020). It has been widely cultivated in south-central and southern China, such as Hunan, Jiangxi, Guangxi, Zhejiang, Fujian, and Hainan provinces (Liu et al., 2018). The area of which C. oleifera plantation approximately reached 4.4 million hectares with an annual output of over 2.6 million tons of seeds consequently oil yields more than 0.65 million tons (Chen et al., 2020).

Drought stress is one of the most crucial environmental factors impairing photosynthesis, thereby limiting plant growth and yield (Muhammad et al., 2020). Accordingly, drought stress is regarded as one of the main threats to food security in the prevailing climate change era (Siddique et al., 2016). Plants have developed strategies to cope with drought stress, to ensure survival under severe drought stress, which includes stomata closure, osmotic adjustment, and enhanced tolerance level (Muhammad et al., 2020). Benefiting from recent advances in high-throughput sequencing, research to identify drought-related genes and miRNA regulatory- networks to understand various drought response mechanisms has made great progress in plants, such as cowpea (Mishra et al., 2021), durum wheat (Liu et al., 2020), and tobacco (Chen et al., 2017). MicroRNAs (miRNAs) are a class of small RNAs that are increasingly being recognized as important modulators of gene expression at the post-transcriptional level. Many miRNAs are involved in drought stress responses, including ABA response, auxin signaling, osmoprotectant, and antioxidant defense, by downregulating the respective target genes encoding regulatory and functional proteins (Ding et al., 2013).

Although C. oleifera has been classified as a drought-tolerant tree species, it is mainly planted on hills and mountains distributed in south China, which has severe difficulties in irrigation. In addition, the southern monsoon climate region has high temperatures and little rainfall in summer, resulting in drought stress becoming a major obstacle to the development of the C. oleifera planting industry.

This study aimed to identify the main candidate genes and miRNAs associated with the different responses between two genotypes of C. oleifera cultivars when they suffered from severe drought stress and to preliminarily reveal the potential molecular mechanisms of miRNA-mRNA regulatory networks that participated in the formation of drought tolerance.

Materials and methods

Plant material and treatments

Two Oil-tea cultivars, namely, C. oleifera cv. ‘Xianglin No. 1’ (XL1) and ‘Hengdong No. 2’ (HD2) were used in this study. The two cultivars were selected by applying a previous study in which twelve C. oleifera cultivars were subjected to severe artificial drought stress until death due to water cutoff. The growth status of the seedlings was observed every 3 days to evaluate the degree of drought stress. After 18 days of continuous water withholding treatment, the seedling survival rate of XL1 was about 80%, while it was zero for HD2, which showed the most obvious difference in survival rate, reflecting that XL1 was more tolerant to drought stress than HD2. The artificial simulated drought stress was conducted by withholding water to an extreme soil water content of 25 percent of the field capacity for 12 days with 3-year-old seedlings, while the control group of seedlings was irrigated daily to normal field capacity, in the mid to late July 2021 in the glasshouse at Hunan Academy of Forestry Sciences which geographical coordinates are 113°01′ east longitude and 28°06′ north latitude. Leaves of two cultivars at 0, 4, 8, and 12 days after drought stress were used for photosynthetic physiological traits measurement and high throughput sequencing library construction. Each cultivar in this experiment has three biological replicates.

Measurement of the photosynthetic physiological traits

The gas exchange indicators include the net photosynthetic rate (P_n, μmol^⋅m^–2s^–1), stomatal conductance (G_s, mmol^⋅m^–2s^–1), intercellular CO₂ concentration (C_i, μmol^⋅mol^–1), and net transpiration rate (T_r, mmol^⋅m^–2s^–1) were measured on the morning during drought stress treatment by using a gas-exchange system LI-6400XT (LI-COR, Lincoln, NE, USA) with the parameter settings were 1,000 μmol photons m^–2⋅s^–1 and 400 CO₂ μmol^⋅mol^–1 (He et al., 2021). Dark-adapted maximum photochemical efficiency of PSII (Fv/Fm) was measured with shaded leaves of each seedling.

RNA libraries construction and sequencing

Twenty-four leaf samples collected from seedlings of the two C. oleifera cultivars that suffered drought stress with three biological replicates were used for RNA library construction and deep sequencing analysis. The mRNA and miRNA library preparations were analyzed on the Illumina Hiseq 2500 platform (Illumina, San Diego, CA, USA) at Allwegene (Beijing, China) with pair-end (2 × 150 bp) and single-end (50 bp) strategies, respectively. Clean data from sequencing were obtained by conducting a series of standardized processing of the raw data, then a certain length range from clean reads was chosen to do all the downstream analysis.

Target prediction and functional annotation

The clean reads obtained from miRNA-seq were mapped to the reference sequence by Bowtie (Langmead et al., 2009) without mismatch to analyze their expression and distribution on the reference. miRBase20.0 was used as a reference to obtain the potential miRNA with modified software mirdeep2 and srna-tools-cli (Friedlander et al., 2012). The available software miREvo (Wen et al., 2012) and mirdeep2 (Friedlander et al., 2012) were integrated to predict novel miRNA. Predicting the target gene of miRNA was performed by psRobot (Wu et al., 2012) for plants. All these analyses were run automatically with default parameters.

Gene Ontology (GO) (Young et al., 2010) and KEGG (Kanehisa et al., 2008) pathway enrichment analysis were conducted on the expressed genes identified from mRNA-seq and target gene candidates of expressed miRNAs identified from miRNA-seq.

Differentially expressed genes and microRNAs analysis

Clean data from mRNA-seq were aligned to the reference genome (Camellia sinensis genome GCF_004153795.1) by TopHat (v2.1.0). HTSeq was used to calculate the number of clean reads aligned to the characterized gene loci, and DESeq2 was then applied to identify the differentially expressed genes (DEGs) with cutoff | log₂ fold change| ≥ 2 and p-value ≤ 0.05. To identify drought-responsive miRNAs, we selected miRNAs with TPM > 1 in at least one of the sample pairs to be analyzed. The criteria to consider a miRNA to be drought-responsive are | log₂ fold change| ≥ 1 and p-value ≤ 0.05.

Co-expression network analysis for construction of modules

We conducted the weighted gene co-expression network analysis (WGCNA) with R package (v1.71) (Langfelder and Horvath, 2008). A total of 2,550 genes that differentially expressed among XL1 and HD2 during different time points of artificial simulation drought stress were selected for further analysis. The gene expression values were imported into WGCNA to construct co-expression modules using the automatic network construction function blockwise modules with default settings, except that the power is 6, TOM Type is unsigned, mergeCutHeight is 0.60, and minModuleSize is 50. The expression profile of merged modules was summarized by the first principal component (module eigengenes, MEs), and the MEs of merged modules were correlated with the physiological traits to find the key modules associated with drought adaptation potential in C. oleifera. The top ten genes with maximum intra-modular connectivity were considered hub genes (Wisniewski et al., 2013) and visualized by Cytoscape (version 3.9.2) (Shannon et al., 2003).

Real-time quantitative RT-PCR

The gene expression levels of four differentially expressed miRNAs including miR408, miR166, and miR398, together with their predicted target genes that were identified from the differentially expressed genes under drought stress among XL1 and HD2, as well as photosynthesis-related genes, plant hormone signal genes, and circadian clock genes that respond to drought stress were determined by qRT-PCR, within use the leaf samples that were the same as those applied to high throughput sequencing according to the method previously reported (He et al., 2021). The primers used for qRT-PCR were determined on the 5′ region of the candidate genes with putative amplicons of 150–250 bp and listed in Table 1.

TABLE 1

Table 1. Information on primers used in the quantitative real-time PCR (qRT-PCR) analysis.

Data analysis

SPSS was used for statistical analyses of the data of the CO₂ exchange, and the analysis of variance (ANOVA) was conducted by using Duncan’s multiple range test, different lowercase letters indicate significant differences at P < 0.05.

Results

Verification of drought tolerance of the two Camellia oleifera cultivars

To verify the different physiological responses to drought stress between the two cultivars, photosynthetic traits were measured in this study. The results showed that drought stress induced a significant reduction of P_n, G_s, C_i, and T_r in both two cultivars (Figure 1). In the comparison of the two cultivars, XL1 kept higher levels of Fv/Fm after 4, 8, and 12 days of drought stress treatment (0.79, 0.77, and 0.71) than those of HD2 (0.74, 0.71, and 0.66), respectively. The leaf relative water content of XL1 after 4, 8, and 12 days of drought stress treatment (0.51, 0.44, and 0.40) also maintained higher levels compared to those of HD2 (0.43, 0.39, and 0.34), respectively.

FIGURE 1

Figure 1. The differences in photosynthetic parameters between XL1 and HD2 during the drought stress treatment. Net photosynthetic rate (P_n), Stomatal conductance (G_s), Intercellular CO₂ concentration (C_i), and Net transpiration rate (T_r) of both XL1 and HD2 leaves were measured by LI-6400XT between 8 and 11 a.m. local time. (A) Mean of the P_n, (B) mean of the G_s, (C) mean of the C_i, (D) mean of the T_r, (E) mean of the Fv/Fm, and (F) mean of the relative water content. Each data point represents the average of three biological replicates. Error bars represent mean ± SE. Different lowercase letters indicate significant differences at P < 0.05.

Transcriptomic changes responded to drought stress in Camellia oleifera

High-throughput Illumina sequencing yielded a total of 883,425,114 clean reads with Q30 above 91.75% from mRNA libraries of leave samples which were mentioned above. Raw data has been uploaded to the Sequence Read Archive as a BioProject (PRJNA875963). These clean reads were mapped into the Camellia sinensis genome (GCF_004153795.1), of which the mapping rate of each sample in this study was ranging from 83.60 to 90.29%, indicating that the reference genome will meet the research needs, and a total of 13,954 expressed unigenes (Both Read count and FKPM ≥ 1) were obtained. Genes with differential expression patterns between sensitive and tolerant cultivars (cutoff fold change ≥ 2 and p-value ≤ 0.05) were defined, and a total of 1,992 differentially expressed genes (DEGs) during the drought stress treatment, including 496 and 1,768 DEGs were detected in XL1 and HD2, respectively (Figure 2). Meanwhile, a total of 1,526 DEGs at different stages during the artificial simulation drought stress treatment among two different cultivars, including 121, 331, 268, and 1,161 DEGs after 0, 4, 8, and 12 days of drought stress treatment, respectively.

FIGURE 2

Figure 2. The venn diagram of differentially expressed genes (DEGs). (A) DEGs response to drought stress in XL1, (B) DEGs response to drought stress in HD2, (C) DEGs between XL1 and HD2 during drought stress.

Functional annotation and classification of differentially expressed genes

According to GO functional enrichment analysis, more than 33 GO biological process terms were significantly enriched, including biosynthetic, organic cyclic compound biosynthetic, oxidation-reduction small molecule metabolic, lipid metabolic process, and photosynthesis (Figure 3). The cellular component categories of DEGs were significantly enriched for the thylakoid, photosynthetic membrane, photosystem, and photosystem II. Significantly enriched GO molecular terms included oxidoreductase activity, calcium ion binding, obsolete coenzyme binding, and antioxidant activity.

FIGURE 3

Figure 3. Gene ontology (GO) analysis for significant DEGs. Graphical representation of numbers of DEGs enriched in biological process, cellular component, and molecular function GO-terms. The symbol “*” represents the significantly enriched GO terms.

The DEGs were associated with various KEGG pathways involved in Photosynthesis, Photosynthesis - antenna proteins, Carbon fixation in photosynthetic organisms, monoterpenoid biosynthesis, circadian rhythm – plant, glyoxylate and dicarboxylate metabolism, glutathione metabolism, porphyrin and chlorophyll metabolism, carbon metabolism and plant hormone signal transduction (Figure 4).

FIGURE 4

Figure 4. Pathway analysis for significant DEGs. Graphical representation of numbers of DEGs enriched in common KEGG pathway.

Differentially expressed genes associated with photosynthesis-related pathways

Genes involving the photosynthesis-related pathways were divided into two groups (Figure 5). The majority of genes involved in photosynthesis and photosynthesis – antenna proteins were enriched in the first group. On the other hand, most of the genes involved in carbon fixation in the photosynthetic organisms pathway were enriched in the second group.

FIGURE 5

Figure 5. The heatmap of DEG involve in photosynthesis-related pathways. The heatmap was generated from the log₂ (FPKM + 1) of three replicates for each sample, and the values of clustering distance were “correlation” and the clustering algorithms was “complete.”

There were 9 DEGs encoded photosystem I reaction center subunit and 9 genes encoded photosystem II reaction center subunit, which down-regulated during the drought stress, meanwhile keeping higher expression levels in XL1 than that in HD2. There were 8 genes involved in photosynthetic electron transportation, including PetC, PetE, PetF, and PetH, whose expression patterns also showed that higher in XL1 than that in HD2 after drought stress.

There were 5 and 9 DEGs encoded proteins constructing LHCI and LHCII, respectively. The expression levels of Lhca3, Lhca4, Lhca5, Lhcb1, Lhcb2, Lhcb3, Lhcb5, and Lhcb7 were up-regulated in XL1 at the first stages of drought stress, while those in HD2 were down-regulated.

Ribulose 1, 5-bisphosphate carboxylase/oxygenase (Rubisco), phosphoglycerate kinase (PGKs), glyceraldehyde 3-phosphate dehydrogenase (GADPHs), and fructose-bisphosphate aldolase (FBAs) played important roles in the CO₂ fixation of the C3 plants. Two genes encoding Rubisco small chain were down-regulated in both XL1 and HD2 during drought stress with relatively higher expression levels in XL1 than that in HD2. Three PGK genes as well as the GADPH gene maintained a relatively stable expression pattern in XL1, whereas these genes were down-regulated in HD2 during drought stress. There were four genes encoding FBAs were identified, among these three genes were down-regulated in the first stages in XL1 and then up-regulated, whereas these genes were up-regulated in HD2.

Phosphoenolpyruvate carboxylase (PEPCs), phosphoenolpyruvate carboxykinase (PEPCKs), alanine transaminase (ALTs), and aspartate aminotransferase (AATs) played important role in the CO₂ fixation in C4-dicarboxylic acid cycle. The PEPC2 gene was up-regulated in XL1 and down-regulated in HD2 after drought stress. There were two genes encoding PEPCKs identified from DEGs, one of those was up-regulated in both XL1 and HD2, and the other was down-regulated in both XL1 and HD2. One ALT1 gene and one of the ALT2 genes were up-regulated in XL1 after drought stress, whereas these two genes were down-regulated in HD2. The rest ALT2 gene maintained a relatively stable expression pattern in XL1, while down-regulated in HD2. The ATT genes maintained a relatively stable expression pattern in XL1, while up-regulated in HD2.

Differentially expressed genes associated with circadian clock system

There were 15 DEGs enriched in the circadian rhythm – plant pathway (Figure 6), among these genes the Cryptochrome (CRY) and Early Flowering 3 (ELF3) encoding proteins participated in the input pathway, while the Late Elongated Hypocotyl (LHY), Timing of CAB expression 1 (TOC1), pseudo-response regulator 5 (PRR5) and GIGANTEA (GI) encoding proteins that constructed core feedback loop in the oscillator.

FIGURE 6

Figure 6. The heatmap of differentially expressed genes (DEGs) involve in the Circadian rhythm - plant pathway. The heatmap was generated from the log₂ (FPKM + 1) of three replicates for each sample, and the values of clustering distance were “correlation” and the clustering algorithms was “complete.”

The gene expression levels of CRY1.1 were up-regulated after drought stress in both XL1 and HD2, while its downstream regulatory gene ELF3 was down-regulated. The LHY gene expression was up-regulated in XL1 and down-regulated in HD2 after drought stress, respectively, while the expression levels of its repressor gene TOC1 were down-regulated in both XL1 and HD2. The expression levels of APRR5 and GI showed no significant differences at the first stages after drought stress in XL1, which showed significantly down-regulated in HD2. The Chalcone Synthase (CHS) acted as an output pathway regulator, there were three CHS genes differentially expressed during drought stress between XL1 and HD2, among these genes down-regulated after 4 days in XL1 and 8 days in HD2, respectively.

Differentially expressed genes associated with phytohormone signal transduction

A total of 32 DEGs involved in several plant hormone signal transduction pathways were identified (Figure 7), including auxin, abscisic acid (ABA), brassinosteroid (BR), gibberellin (GA), salicylic acid (SA), ethylene (ET), jasmonic acid (JA) and cytokinin (CK), respectively. There were 8 genes involved in the auxin signaling pathway, among these genes 5 genes encoding auxin-responsive protein IAAs and 3 genes encoding auxin-responsive protein SAURs, respectively. The gene expression levels of three IAAs genes showed up-regulated at the first stages after drought stress and two IAAs gene were down-regulated in XL1, while the expression levels of all the five IAAs genes showed down-regulated in HD2. The ABA signaling pathway consisted of the abscisic acid receptor (PYR/PYL), protein phosphatase 2C (PP2C), sucrose non-fermenting 1-related protein kinase 2 (SnRK2), and their downstream regulatory ABA-responsive element binding factor (ABF). The expression level of PYL6 gene which encodes an abscisic acid receptor was up-regulated in both cultivars after drought stress, while the expression levels in XL1 were higher than those in HD2 at each time point under drought stress, respectively. There were four genes encoding PP2Cs family protein, including PP2C16, PP2C24, PP2C51, and PP2C75, among these genes PP2C24, PP2C51, and PP2C75 gene expressions were inhibited by drought stress in XL1, while PP2C16, PP2C24, and PP2C51 gene expressions were induced by drought stress in HD2. The expression levels of SnRK2.8 gene were significantly down-regulated at first in both XL1 and HD2, and then slightly recovered. The bZIP8 and ABI5 genes as downstream regulatory genes in the ABA signaling transduction pathway, their expression levels increased to a certain degree under drought stress, which indicated that the ABA signaling transduction pathway played a crucial role in the molecular response to drought stress. Another plant hormone signal transduction pathway participated in drought adaptation regulation involved with BR, the BSK4 and BSK8 genes as the BR signal receptor kinase encoding genes, their expression levels increased to a certain degree under drought stress, which indicated that the BR signal transduction pathway might be activated by drought stress. The GSK2 and BZR1 genes were up-regulated and down-regulated after drought stress in both cultivars, respectively, consequently, the downstream regulatory CYCD3-1 gene was down-regulated after drought stress. In the SA pathway, the NPR5 and TGA2 genes were both up-regulated during drought stress in XL1 and HD2, while the downstream regulatory PR-1 gene was down-regulated in XL1 and up-regulated in HD2 during drought stress, respectively. The key regulatory genes DELLA from GA pathway and EIL from ET pathway were both induced during drought stress in XL1 and HD2, indicating that these pathways might participate in downstream drought-responsive processes. The above results suggested that the plant hormone signal transduction pathway played important role in drought adaptation in C. oleifera.

FIGURE 7

Figure 7. The heatmap of differentially expressed genes (DEGs) involve in the Plant hormone signal transduction pathway. The heatmap was generated from the log₂ (FPKM + 1) of three replicates for each sample, and the values of clustering distance were “correlation” and the clustering algorithms was “complete.”

Drought stress responsive transcription factors

In this study, a total of 150 transcription factors (TFs) were differentially expressed during the drought stress treatment. The most enriched TF families that responded to drought stress were bHLH, bZIP, C2H2, ERF, GATA, GRAS, HD-ZIP, MYB, MYB related, NAC, TCP, and WRKY. The MYB and MYB-related family was the largest group (12.00%), followed by the bHLH family (8.00%), WRKY family (8.00%), and NAC family (7.33%). There were 6 and 10 genes differentially expressed during the whole drought stress treatment in XL1 and HD2, respectively (Figure 8). The heatmap of differentially expressed TFs showed that there were 8 TFs significantly up-regulated including two from the HSF family (HSF8 and HSF24), and one from each of MYB related, bHLH, NAC, bZIP, C3H and ERF families, while the rest of differentially expressed TFs significantly down-regulated including TFs that mainly from the bZIP, NAC, LBD, bHLH, and SBP families.

FIGURE 8

Figure 8. The heatmap of differentially expressed transcription factors (TFs) during the drought stress treatment. The heatmap was generated from the log2 (FPKM + 1) of three replicates for each sample, and the values of clustering distance were “correlation” and the clustering algorithms was “complete.”

Weighted gene co-expression network analysis of differentially expressed genes

The WGCNA result (Figure 9) indicated that DEGs could be grouped into 6 modules, which contained 1 327, 574, 285, 210, 77, and 77 in turquoise, brown, blue, black, magenta, and pink, respectively. Correlation analysis that conducted between physiological traits and module eigengenes (MEs) indicated that genes of brown module were significantly positive correlated to the photosynthetic physiological traits such as P_n (r² = 0.73, p < 0.05), G_s (r² = 0.80, p < 0.05) and Fv/Fm (r² = 0.74, p < 0.05), as well as RWC (r² = 0.85, p < 0.05). In contrast, genes clustered in turquoise module negatively correlated to P_n (r² = −0.75, p < 0.05), G_s (r² = −0.61, p < 0.05) and Fv/Fm (r² = −0.80, p < 0.05), as well as RWC (r² = −0.71, p < 0.05).

FIGURE 9

Figure 9. Network analysis dendrogram showing modules identified by WGCNA. (A) Hierarchical cluster tree showing 9 co-expression modules (including gray module) identified by WGCNAs, each branch in the tree represents an individual gene; (B) Correlation between module eigengenes and photosynthetic parameters, each cell contains the corresponding correlation and p-value.

Brown and turquoise modules oppositely and positively correlated with photosynthetic physiological traits were selected to construct Protein-Protein Interaction (PPI) networks by enrichment analysis through the STRING database with the predicted proteins encoded by genes from the two modules (Figure 10). Functional enrichment in the PPI network showed that genes in brown and turquoise modules enrich in photosynthetic electron transport in photosystem I, negative regulation of response to salt stress, negative regulation of response to water deprivation, photosystem II assembly and gluconeogenesis Biological Process; chlorophyll-binding, rRNA binding, isomerase activity, calcium ion binding and oxidoreductase activity Molecular Function; photosystem II oxygen-evolving complex, photosystem II, photosystem and chloroplast thylakoid lumen Cellular Component. The KEGG pathways enriched with genes in brown and turquoise modules contained photosynthesis – antenna proteins, photosynthesis, thiamine metabolism, glyoxylate and dicarboxylate metabolism, and porphyrin and chlorophyll metabolism. These results suggest that photosynthesis-related genes were significantly affected during the drought stress, their differences between drought-tolerant and sensitive cultivars might play a crucial role in regulating drought adaptation in C. oleifera.

FIGURE 10

Figure 10. Protein-Protein Interaction network analysis of drought-responsive modules. PPI networks of the positively correlated brown module and the negatively correlated turquoise module were visualized using the Cytoscape software platform. The correlation networks of the top 10 genes and their directly regulatory genes from brown and turquoise modules were clustered into two circles. The color depth of 10 hub genes represents the strength of connectivity.

To screen for candidate hub genes, the top ten connected genes in the selected modules were ranked by Maximal Clique Centrality (MCC) method with Cytoscape (Figure 11), resulting in 3 genes from the brown module and 7 genes from the turquoise module were screened out, respectively. In the related brown module, the hub genes encoded Thylakoid lumenal 16.5 kDa protein, F-type h + -transporting ATPase subunit delta, and Phosphoribulokinase. In the related turquoise module, the hub genes encoded Thioredoxin-dependent peroxiredoxin, Thylakoid lumen 18.3 kDa family protein, Sedoheptulose-1,7-bisphosphatase, Curvature Thylakoid 1c, Thylakoid membrane phosphoprotein 14 kDa, Photosynthetic ndh subunit of subcomplex b 5 and an Uncharacterized protein. This indicates that these genes may play a crucial role in the regulatory network of drought stress adaptation.

FIGURE 11

Figure 11. Expression levels of top ten hub genes of the PPI network. Heatmap showing the expression profiles of the top ten hub genes, the color depth of each cell represents the gene expression level by log₂(FPKM + 1), and the values of clustering distance were “correlation” and the clustering algorithms was “complete.”

Identification of miRNAs in response to drought stress in Camellia oleifera

High-throughput Illumina sequencing yielded a total of 52,492,335 unique reads in the miRNA libraries, and sequences with lengths ranging from 18 to 35 nt were selected for the following analysis. Mapped small RNA tags were used to look for known miRNAs against miRase and 3,051 known miRNAs from 64 miRNA families were identified. In addition, de novo predictions have been performed that identified another 241 miRNA candidates.

To explore the differential expression of miRNAs related to drought stress in XL1 and HD2, 83 expressed miRNAs (readcount ≥ 1 in each sample), including 55 known miRNAs and 28 novel miRNAs, were conducted differential expression analysis by DESeq2 package. Comparisons with control groups revealed 7 and 5 differential expressed miRNAs (DEMs) in XL1 and HD2, respectively (cutoff fold change ≥ 1 and p-value ≤ 0.05). Among these DEMs, 5 were identified as known miRNA mainly from miR398 and miR408 families, while 2 were predicted as novel miRNA candidates. Known and novel miRNAs with differential expression patterns between sensitive and tolerant cultivars were defined, a total of 14 miRNAs during the artificial simulation of drought stress treatment, including 0, 6, 6, and 5 miRNAs after 0, 4, 8, 12 days of drought stress treatment, respectively, were identified as DEMs between two cultivars. Among these DEMs, 11 were identified as known miRNA from miR159, miR166, miR395, miR398, miR408, miR482, miR4995, and miR5368 families, while 4 were predicted as novel miRNAs.

Target prediction and construction of a regulatory and interaction network

Differentially expressed miRNAs target gene prediction and GO annotation analysis was performed to characterize the regulatory roles of miRNAs in response to drought stress, and a total of 78 target DEGs were found. For these target DEGs, among the GO terms identified in the biological process, eight significantly enriched categories were “generation of precursor metabolites and energy,” “macromolecule localization,” “cellular component organization,” “cellular component organization or biogenesis,” “organic substance transport,” “cellular protein localization,” “cellular macromolecule localization” and “response to stress.”

To establish the regulatory network of miRNA-mRNA involved in the response to drought stress, the predicted target DEGs of DEMs were analyzed (Figure 12). The correlation analysis of expression levels was conducted between DEMs and the predicted target DEGs, 33 miRNA-mRNA pairs were screened out, involving 9 DEMs and 32 DEGs. Among these pairs, 12 pairs showed antagonistic regulatory patterns due to their significantly negative correlation-ship with expression levels (p < 0.05) (Table 2).

FIGURE 12

Figure 12. miRNA-mRNA regulatory network involved in the response to drought stress. The network of differentially expressed miRNAs and the predicted target genes was shown, and the negative and positive correlations between DEMs and their target were reflected by the solid line with a stop symbol and dotted line with an arrow, respectively.

TABLE 2

Table 2. miRNA-mRNA regulatory network response to drought stress in Camellia oleifera.

Verification of differentially expressed genes and differentially expressed miRNAs analysis results by quantitative real-time PCR

To validate the reliability of the miRNA-seq and mRNA-seq data, 8 DEGs and 4DEMs that showed significant differences in expression were selected: RbcS and FBA from carbon fixation in photosynthetic organisms pathway (Figures 13A,B), PLY6 and DELLA from plant hormone signal transduction pathway (Figures 13C,D), LHY and TOC1 circadian rhythm - plant pathway (Figures 13E,F), miR398, miR408-3p, miR166 (Figures 13G–I), and their predicted miRNA candidates (Figures 13J–L). We analyzed the variance of expression of DEGs and DEMs in samples during the whole drought stress treatment using RT-PCR. The mRNA-seq and RT-PCR data were very closely correlated, and there was high consistency in the up- and down-regulated expression of DEGs and DEMs.

FIGURE 13

Figure 13. The comparison of the expression levels of 6 DEGs, 3 DEMs and their target genes identified response to drought stress among mRNA-Seq/miRNA-seq (line) and qRT-PCR analysis (bar). (A) Rbcs gene expression levels during the drought stress by qRT-PCR and mRNA-seq quantifications. (B) FBA gene expression levels during the drought stress by qRT-PCR and mRNA-seq quantifications. (C) PYL6 gene expression levels during the drought stress by qRT-PCR and mRNA-seq quantifications. (D) DELLA gene expression levels during the drought stress by qRT-PCR and mRNA-seq quantifications. (E) LHY gene expression levels during the drought stress by qRT-PCR and mRNA-seq quantifications. (F) TOC1 gene expression levels during the drought stress by qRT-PCR and mRNA-seq quantifications. (G) miR398 gene expression levels during the drought stress by qRT-PCR and mRNA-seq quantifications. (H) miR408-3p gene expression levels during the drought stress by qRT-PCR and mRNA-seq quantifications. (I) miR166 gene expression levels during the drought stress by qRT-PCR and mRNA-seq quantifications. (J) DMR6 gene expression levels during the drought stress by qRT-PCR and mRNA-seq quantifications. (K) EDR2 gene expression levels during the drought stress by qRT-PCR and mRNA-seq quantifications. (L) LOC114322616 gene expression levels during the drought stress by qRT-PCR and mRNA-seq quantifications.

Discussion

Photosynthetic physiology traits response to drought stress in Camellia oleifera

Photosynthesis-related pathways consisted of photosynthetic pigments and photosystems, electron transport chain, and carbon dioxide reduction pathways. Damages to any of these components reduce the overall synthetic capacity of plants (Siddique et al., 2016). Plants perceive drought through chemical and hydraulic signaling and trigger chemical signaling from the root toward the shoot to initiate necessary molecular responses (Mukarram et al., 2021). These responses induce morpho-physiological changes, including stomatal behavior, and prepare the plant for the upcoming drought (Janiak et al., 2016). In general, dynamic regulation of the stomatal aperture prevent excessive transpirational water loss and minimizes the drought impact (Lawson and Matthews, 2020). Stomatal conductance provides an easy, integrated measure of the physiological response of plants to drought stress (Gilbert and Medina, 2016). The stomatal closure is triggered when plants suffered water deficiency, consequently, plants prioritize survival over growth and development. In this study, the decrease of G_s soon after drought stress in XL1 and HD2 reflected that both cultivars initiated the active closure mechanism to prevent water loss for further damage. Plants sacrifice photosynthetic capacities for survival when they undergo severe drought stress in the way of active stomatal closure that leads to less uptake CO₂ abilities. In the comparison of the two cultivars, the G_s remain at higher levels in XL1 than that in HD2, while the P_n and T_r similarly reduced within the deepening drought stress with consistently higher levels in XL1 than that in HD2, indicating that XL1 remains higher photosynthetic capacities than HD2. Dark-adapted maximum photochemical efficiency of PSII (Fv/Fm) is a useful indicator of past stress in leaves and acclimation. Water deficiency limits the efficacy of the photosynthetic apparatus, causes substantial damage to the thylakoid membrane, and reduces the Chl contents of leaves, thus affecting the photoinhibition/photochemical intensity, resulting in reduced Fv/Fm (Muhammad et al., 2020). Fv/Fm depression is sustained and may decrease photosynthetic rates and thus also represent damage (Demmig-Adams et al., 2012). In this study, the Fv/Fm ratio significantly decreased in both cultivars after drought stress with a relatively higher level in comparison to XL1 and HD2, which indicated that HD2 suffered more severe damages than XL1 within the photosystem and electron transport chain.

Differentially expressed transcription factors under drought stress in Camellia oleifera

Transcription factors (TFs), as the binding to the cis-regulatory elements in promotors to mediate the expression of massive genes, have fundamental importance in critical aspects of plant function in ensuring that plant growth and development match fluctuations in environmental conditions (Hrmova and Hussain, 2021). Notably, fully sequenced genomes of several plants showed that around 10% of all identified genes encode TFs (Baillo et al., 2019). In this study, a large number of DEGs encoding members of TF families including bHLH, bZIP, MYB, NAC, and WRKY, were identified according to their significantly responsive expression changes during drought stress. MYB TFs as one of the largest protein families in plants have been found to participate in diverse processes such as plant development and responses to biotic and abiotic stresses (Guo et al., 2022). In this study, we have found a drought-responsive differentially expressed gene encoding Reveille1 (RVE1), a Myb-like, clock-regulated TF that is homologous to the central clock genes Circadian Clock Associated 1 (CCA1) and Late Elongated Hypocotyl (LHY). A previous study revealed that RVE1 acts as a clock output instead of a central oscillator, due to the disordered expression of RVE1 does not affect circadian rhythmicity (Rawat et al., 2009). The Phytochrome Interacting Factor (PIF) is a small subset of basic helix-loop-helix TFs that function as a cellular signaling hub, integrating internal and external stimuli to coordinate the regulation of the transcriptional network that drives multiple morphogenic responses (Leivar and Quail, 2011). The PIF4, PIF5, and PIF7 have been shown that act as clock output controlled by circadian clock genes that manipulated plant growth and development (Cheng et al., 2021; Favero et al., 2021). It is reported that PIF7 negatively regulates DREB1 expression under circadian control in Arabidopsis, consequently repressed DREB1 expression may be important for avoiding plant growth retardation under unstressed conditions. In this study, a PIF7 gene significantly decreased in HD2 under drought stress, while it has no significant differences in XL1. These indicate that the circadian clock system might involve in the differences in drought stress response.

miRNA-mRNA network under drought stress in Camellia oleifera

Plants have developed the function of using miRNA to regulate gene transcription in response to environmental factors and stress conditions such as temperature, light, water, and nutrient deficiency (Kryovrysanaki et al., 2022). Among the screened differentially expressed miRNAs in this study, miR398 played an important role in the regulatory network that responds to drought stress in C. oleifera. Several studies have been made for the identification of drought-responsive miRNAs, and members from the miR398 family have been proven to be one of the most important stress-responsive miRNAs, which play important roles in oxidative stress, drought stress, salt stress, ABA signal, copper and phosphate deficiency, as well as biotic stress (Bakhshi et al., 2014). Quantitative expression analysis of miRNAs during drought stress in Camellia sinensis carried out with the most famous oolong tea cultivar ‘Tieguanyin’ by former researchers revealed that members from miR166, miR398, and miR408 families were possibly associated with the C. sinensis drought stress (Guo et al., 2017), which provided evidence that supported the credibility of our results. Interestingly, the expression levels of miR166 and miR398 were down-regulated and up-regulated during the different phases of drought stress, respectively, which were exactly the opposite of those in C. oleifera from our results. Commonly, the network of miRNAs could be the same in different plant species, or it may differ between species (Akdogan et al., 2016). The consistent results were obtained from research conducted in Echinacea purpurea L., which showed that miR398 expression was reduced in comparison with the control in different field capacities of water (Moradi and Khalili, 2019). Furthermore, a previous study illustrated that miR398 targets Cu/Zn-superoxide dismutase (CSD1 and CSD2) which were found to be significantly induced in the leaf of wheat cultivars to provide drought tolerance (Akdogan et al., 2016). In this study, the expression level of miR398b negatively correlated to those of gene-LOC114259217, gene-LOC114291719, gene-LOC114268758, and gene-LOC114287709, which encoding a DOWNY MILDEW RESISTANCE 6-like protein, a phosphoinositide phosphatase SAC1-like protein, a mitochondrial import receptor subunit TOM40-1-like protein and a glutamate/aspartate-prephenate aminotransferase-like protein. The significant decreases of miR398b were observed and consequently led to increasing the expression levels of the predicted target genes mentioned above. Besides, the similar expression between other predicted target genes and miRNAs indicated the inconsistent regulation of miRNA, these might be mostly due to the expression of unknown TFs and activation of other metabolic networks in the drought stress resistance network (Moradi and Khalili, 2019). In this study, these unexpected regulatory networks are mainly involved in gene-LOC114257452, gene-LOC114320625, and gene-LOC114291874, which encode THO complex subunit 4A, Chloroplast Unusual Positioning 1 protein, 60S ribosomal protein L44-like.

miR408 is a multistress-responsive miRNA that observed to be differentially regulated in plants such as Arabidopsis (Song et al., 2017), Triticum (Feng et al., 2013), Oryza (Balyan et al., 2022), Chickpea (Hajyzadeh et al., 2015), and Vigna unguiculata (Mishra et al., 2021). Previous studies showed that overexpression of miR408, by down-regulation of the direct target genes, obviously promotes photosynthetic activity, biomass, and seed yield (Song et al., 2017; Pan et al., 2018), as well as tolerance to drought stress (Hang et al., 2021). However, another study of overexpression mutant revealed that higher miR408 expression led to enhanced sensitivity to drought stress (Hajyzadeh et al., 2015; Ma et al., 2015). This opposite role of miR408 in plant drought-responsive regulatory networks might indicate a species-specific drought-responsive characteristic for miR408 (Mutum et al., 2013). The miR408b-3p identified in our study showed that it was down-regulated after drought stress, and the abundance in XL1 was consistently higher than that in HD2. The miR408-3p induced its predicted target gene which encodes protein Enhanced Disease Resistance 2 by negatively regulating of decreasing itself, might further improve drought tolerance via crosstalk between different stress-responsive pathways.

Data availability statement

The datasets presented in this study can be found in online repositories. The name of the repository and accession number can be found below: NCBI Sequence Read Archive; PRJNA875963.

Author contributions

ZH and CL: software and writing—review and editing. ZZ: validation and data curation. RW and ZH: investigation. ZH: writing—original draft preparation. YC: supervision, project administration, and funding acquisition. All authors read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Key R&D Program, China (Grant No. 2018YFD1000603) and the Natural Science Foundation of Hunan Province, China (Grant No. 2022JJ30325).

Acknowledgments

We thank Allwegene for the efficiently high throughput sequencing project.

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.

References

Akdogan, G., Tufekci, E. D., Uranbey, S., and Unver, T. (2016). miRNA-based drought regulation in wheat. Funct. Integr. Genomics 16, 221–233.

Google Scholar

Baillo, E. H., Kimotho, R. N., Zhang, Z., and Xu, P. (2019). Transcription Factors Associated with Abiotic and Biotic Stress Tolerance and Their Potential for Crops Improvement. Genes 10:771.

Google Scholar

Bakhshi, B., Fard, E. M., and Hosseini, G. (2014). Evaluation of miR398 differential expression in rice under drought stress condition. Bull. Georg. Natl. Acad. Sci. 10:771.

Google Scholar

Balyan, S., Kansal, S., Jajo, R., Behere, P. R., Chatterjee, R., and Raghuvanshi, S. (2022). Delineating the Tissue-Mediated Drought Stress Governed Tuning of Conserved miR408 and Its Targets in Rice. Preprints 2022:2022020113. doi: 10.20944/preprints202202.0113.v1

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Q., Li, M., Zhang, Z., Tie, W., Chen, X., Jin, L., et al. (2017). Integrated mRNA and microRNA analysis identifies genes and small miRNA molecules associated with transcriptional and post-transcriptional-level responses to both drought stress and re-watering treatment in tobacco. BMC Genomics 18:62. doi: 10.1186/s12864-016-3372-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Deng, S., Chen, L., Ma, L., He, H., Wang, X., et al. (2020). New Theory on the Development of Camellia Oil. J. Nanjing For. Univ. 44, 657–667.

Google Scholar

Cheng, M. C., Kathare, P. K., Paik, I., and Huq, E. (2021). Phytochrome Signaling Networks. Annu. Rev. Plant Biol. 72, 217–244.

Google Scholar

Demmig-Adams, B., Cohu, C. M., Muller, O., and Adams, W. W. (2012). Modulation of photosynthetic energy conversion efficiency in nature: From seconds to seasons. Photosynth. Res. 113, 75–88. doi: 10.1007/s11120-012-9761-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, Y., Tao, Y., and Zhu, C. (2013). Emerging roles of microRNAs in the mediation of drought stress response in plants. J. Exp. Bot. 64, 3077–3086.

Google Scholar

Favero, D. S., Lambolez, A., and Sugimoto, K. (2021). Molecular pathways regulating elongation of aerial plant organs: a focus on light, the circadian clock, and temperature. Plant J. 105, 392–420. doi: 10.1111/tpj.14996

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, H., Zhang, Q., Wang, Q., Wang, X., Liu, J., Li, M., et al. (2013). Target of tae-miR408, a chemocyanin-like protein gene (TaCLP1), plays positive roles in wheat response to high-salinity, heavy cupric stress and stripe rust. Plant Mol. Biol. 83, 433–443. doi: 10.1007/s11103-013-0101-9

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilbert, M. E., and Medina, V. (2016). Drought Adaptation Mechanisms Should Guide Experimental Design. Trends Plant Sci. 21, 639–647.

Google Scholar

Guo, L., Shen, J., Zhang, C., Guo, Q., Liang, H., Hou, X., et al. (2022). Characterization and bioinformatics analysis of ptc-miR396g-5p in response to drought stress of Paeonia ostii. Non Coding RNA Res. 7, 150–158. doi: 10.1016/j.ncrna.2022.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, Y., Zhao, S., Zhu, C., Chang, X., Yue, C., Wang, Z., et al. (2017). Identification of drought-responsive miRNAs and physiological characterization of tea plant (Camellia sinensis L.) under drought stress. BMC Plant Biol. 17:211. doi: 10.1186/s12870-017-1172-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Hajyzadeh, M., Turktas, M., Khawar, K. M., and Unver, T. (2015). miR408 overexpression causes increased drought tolerance in chickpea. Gene 555, 186–193. doi: 10.1016/j.gene.2014.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Hang, N., Shi, T., Liu, Y., Ye, W., Taier, G., Sun, Y., et al. (2021). Overexpression of Os-microRNA408 enhances drought tolerance in perennial ryegrass. Physiol. Plant 172, 733–747. doi: 10.1111/ppl.13276

PubMed Abstract | CrossRef Full Text | Google Scholar

He, Z., Liu, C., Wang, X., Wang, R., Chen, Y., and Tian, Y. (2020). Assessment of genetic diversity in Camellia oleifera Abel. accessions using morphological traits and simple sequence repeat (SSR) markers. Breed. Sci. 70, 586–593. doi: 10.1270/jsbbs.20066

PubMed Abstract | CrossRef Full Text | Google Scholar

He, Z., Liu, C., Wang, X., Wang, R., Tian, Y., and Chen, Y. (2021). Leaf Transcriptome and Weight Gene Co-expression Network Analysis Uncovers Genes Associated with Photosynthetic Efficiency in Camellia oleifera. Biochem. Genet. 59, 398–421. doi: 10.1007/s10528-020-09995-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Hrmova, M., and Hussain, S. S. (2021). Plant Transcription Factors Involved in Drought and Associated Stresses. Int. J. Mol. Sci. 22:5662.

Google Scholar

Janiak, A., Kwasniewski, M., and Szarejko, I. (2016). Gene expression regulation in roots under drought. J. Exp. Bot. 67, 1003–1014.

Google Scholar

Kanehisa, M., Araki, M., Goto, S., Hattori, M., Hirakawa, M., Itoh, M., et al. (2008). KEGG for linking genomes to life and the environment. Nucleic Acids Res. 36, D480–D484.

Google Scholar

Kryovrysanaki, N., James, A., Bardani, E., and Kalantidis, K. (2022). RNA silencing pathways in plant development and defense. Int. J. Dev. Biol. 66:163.

Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: An R package for weighted correlation network analysis. BMC Bioinform. 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Langmead, B., Trapnell, C., Pop, M., and Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10:R25. doi: 10.1186/gb-2009-10-3-r25

PubMed Abstract | CrossRef Full Text | Google Scholar

Lawson, T., and Matthews, J. (2020). Guard Cell Metabolism and Stomatal Function. Annu. Rev. Plant Biol. 71, 273–302.

Google Scholar

Leivar, P., and Quail, P. H. (2011). PIFs: Pivotal components in a cellular signaling hub. Trends Plant Sci. 16, 19–28. doi: 10.1016/j.tplants.2010.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, C., Chen, L., Tang, W., Peng, S., Li, M., Deng, N., et al. (2018). Predicting Potential Distribution and Evaluating Suitable Soil Condition of Oil Tea Camellia in China. Forests 9:487.

Google Scholar

Liu, H., Able, A. J., and Able, J. A. (2020). Integrated Analysis of Small RNA, Transcriptome, and Degradome Sequencing Reveals the Water-Deficit and Heat Stress Response Network in Durum Wheat. Int. J. Mol. Sci. 21:6017. doi: 10.3390/ijms21176017

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, C., Burd, S., and Lers, A. (2015). miR408 is involved in abiotic stress responses in Arabidopsis. Plant J. 84, 169–187.

Google Scholar

Mishra, S., Sahu, G., and Shaw, B. P. (2021). Integrative small RNA and transcriptome analysis provides insight into key role of miR408 towards drought tolerance response in cowpea. Plant Cell Rep. 41, 75–94. doi: 10.1007/s00299-021-02783-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Moradi, K., and Khalili, F. (2019). Effect of drought stress on expression of HSP70 protein and miR398 in Echinacea purpurea L. J. Med. Herbs 10, 185–194.

Google Scholar

Muhammad, I., Shalmani, A., Ali, M., Yang, Q., Ahmad, H., and Li, F. (2020). Mechanisms regulating the dynamics of photosynthesis under abiotic stresses. Front. Plant Sci. 11:615942. doi: 10.3389/fpls.2020.615942

PubMed Abstract | CrossRef Full Text | Google Scholar

Mukarram, M., Choudhary, S., Kurjak, D., Petek, A., and Khan, M. M. A. (2021). Drought: Sensing, signalling, effects and tolerance in higher plants. Physiol. Plant. 172, 1291–1300.

Google Scholar

Mutum, R. D., Balyan, S. C., Kansal, S., Agarwal, P., Kumar, S., Kumar, M., et al. (2013). Evolution of variety-specific regulatory schema for expression of osa-miR408 in indica rice varieties under drought stress. FEBS J. 280, 1717–1730. doi: 10.1111/febs.12186

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, J., Huang, D., Guo, Z., Kuang, Z., Zhang, H., Xie, X., et al. (2018). Overexpression of microRNA408 enhances photosynthesis, growth, and seed yield in diverse plants. J. Integr. Plant Biol. 60, 323–340. doi: 10.1111/jipb.12634

PubMed Abstract | CrossRef Full Text | Google Scholar

Rawat, R., Schwartz, J., Jones, M. A., Sairanen, I., Cheng, Y., Andersson, C. R., et al. (2009). REVEILLE1, a Myb-like transcription factor, integrates the circadian clock and auxin pathways. Proc. Natl. Acad. Sci. U.S.A. 106, 16883–16888. doi: 10.1073/pnas.0813035106

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, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Siddique, Z., Jan, S., Imadi, S. R., Gul, A., and Ahmad, P. (2016). “Drought stress and photosynthesis in plants,”. Water Stress and Crop Plants: a Sustainable Approach. (ed) P. Ahmad (Hoboken, NJ: John Wiley & Sons) 2, 1–11.

Google Scholar

Song, Z., Zhang, L., Wang, Y., Li, H., Li, S., Zhao, H., et al. (2017). Constitutive Expression of miR408 Improves Biomass and Seed Yield in Arabidopsis. Front. Plant Sci. 8:2114. doi: 10.3389/fpls.2017.02114

PubMed Abstract | CrossRef Full Text | Google Scholar

Wen, M., Shen, Y., Shi, S., and Tang, T. (2012). miREvo: An integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinform. 13:140. doi: 10.1186/1471-2105-13-140

PubMed Abstract | CrossRef Full Text | Google Scholar

Wisniewski, N., Cadeiras, M., Bondar, G., Cheng, R. K., Shahzad, K., Onat, D., et al. (2013). Weighted Gene Coexpression Network Analysis (WGCNA) Modeling of Multiorgan Dysfunction Syndrome after Mechanical Circulatory Support Therapy. J. Heart Lung Transplant. 32:S223.

Google Scholar

Wu, H. J., Ma, Y. K., Chen, T., Wang, M., and Wang, X. J. (2012). PsRobot: A web-based plant small RNA meta-analysis toolbox. Nucleic Acids Res. 40, W22–W28. doi: 10.1093/nar/gks554

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, M. D., Wakefield, M. J., Smyth, G. K., and Oshlack, A. (2010). Gene ontology analysis for RNA-seq: Accounting for selection bias. Genome Biol. 11:R14.

Google Scholar