Transcriptomics analysis of field-droughted pear (Pyrus spp.) reveals potential drought stress genes and metabolic pathways

Plant growth conditions and drought treatment

Field drought experiments were performed for three continuous years in a pear germplasm nursery at the Institute of Fruit, Shanxi Academy of Agricultural Science, beginning on 21 October 2015 and concluding on 21 October 2018. The pear nursery is located in a semi-arid area of Taigu, Shanxi Province, China (37°26′N, 37°26′E) with an altitude of 750 m and managed according to common cultural practices in the region. In this region, the annual average temperature is 9.8 °C with an annual accumulated temperature above 10 °C (AT10) of 3529 °C. The annual hours of sunshine range from 2,500 h to 2,600 h with an average frost-free period of 149 days. The annual rainfall is 450 mm, and the annual accumulative evaporation is 1,800 mm, which is approximately four times higher than the average total rainfall.

The pear cultivar YuluXiangli (Pyrus spp) was used in the experiment. YuluXiangli was derived from a cross between Pyrus bretschneiderie and Pyrus sinkiangensis, and is resistant to drought. The irrigation (control) and field-drought treatments were assigned via a randomized block design with three replicates, where the fields were divided into six plots with 10 healthy and uniform 15-year-old pear trees per plot. Field drought plots were exposed to rainfall without additional irrigation, whereas control plots were irrigated in November, May, and July, each of which received 728.5 tons water/acre. The maximum water holding capacity was 30% in field-drought treatment (severe drought) and 75% to 80% in control with irrigation. Fertilization and pest controls were consistent among the field-drought and control plots.

In total, 100 young leaves, branches, and young fruits, including 10 from each tree, were independently harvested on 5 May 2018, and were swiftly placed in liquid nitrogen and stored at −80 °C for RNA extraction.

At maturity, 10 fruits, each from a single tree, were independently harvested to determine fruit soluble solids content with a handheld PAL-1 digital display sugar meter (Atago, Tokyo, Japan) and single fruit weight.

Total RNA extraction, library preparation and sequencing

Total RNA was extracted from young leaves, branches, and young fruits for each treatment using RNApreo Pure Plant Kit (Tiangen, Beijing, China) in accordance with the manufacturer’s instructions. RNA purity and integrity were determined by Agilent 2100 Bioanalyzer (Agilent Technology, USA) according to the manufacturer’s instructions. The qualified RNA with an RNA integrity number (RIN) of ≥7 and an 28S/18S ribosomal RNA ratio of ≥0.7 was applied to construct 10 cDNA libraries (5 repeats for drought and irrigation, respectively). Equal amounts of RNA from young leaves, branches, and young fruits for each treatment were mixed, and then were diluted to 1 ng/µL for library construction. Briefly, RNA was enriched by magnetic beads containing poly-T oligos and fragmented first to 200–300 bp in length by ion interruption, and reversed transcribed to the first strand of cDNA by 6-bp random primers. Then, the first strand of cDNA was used as a template to synthesize the second strand of cDNA. Library fragments were enriched by PCR amplification to select the fragment size of 300–400 bp. Equal amounts of libraries with different index sequences were pooled prior to sequencing and diluted to 2 nM for paired-end sequencing on the Illumina HiSeq 2500 platform. All raw reads were deposited in the NCBI repository with Bioproject: PRJNA655255 under the accession numbers of SRR12424088 –SRR12424107.

Read mapping and transcript profiling

The adapter and low-quality sequences were removed from the raw RNA-seq reads to generate high-quality clean reads that were aligned to the pear genome reference GCF_000315295.1_Pbr_v1.0_(https://ftp.ncbi.nlm.nih.gov/genomes/all/000/315/295/ GCF_000315295.1_Pbr_v1.0/) with HISAT2 (http://ccb.jhu.edu/software/hisat2/index.shtml). Following the alignments, the raw counts for each pear gene were normalized as fragments per kilobase of transcript per million mapped reads (FPKM) (Trapnell et al., 2010). Principal component analysis (PCA) was performed to compare the log2-transformed FPKM values of the expressed gene profiles among tissue-type and stages using the prcomp function in the R program (https://www.r-project.org/). The hierarchical clustering of samples was performed using Pheatmap in R. Read coverage over gene body was calculated by RSeQC (Wang, Wang & Li, 2012), and the corresponding plot figure was made by using ggplot2 with R script.

Identification of differentially expressed genes (DEGs)

DEGs among tissue-types at different stages were located using the statistical package DEGseq with the MA-plot-based method (Wang et al., 2010) in R version 3.0.3, where genes were considered differentially expressed if —log2FoldChange—>1, and an adjusted p value using Benjamini–Hochberg procedure (Benjamini & Hochberg, 1995) (false discovery rate (FDR)) was <0.05.

Gene annotation (GO) and functional enrichment analysis

The GO enrichment analysis for biological processes, molecular functions, and cellular components was performed using TopGo (Alexa & Rahnenfuhrer, 2016) with P value <0.05. Pathway enrichment analysis was implemented on all DEGs in the Kyoto Encyclopedia of Genes and Genome (KEGG) platform (http://www.genome.jp/kegg/) (Kanehisa et al., 2008). An adjusted P value <0.05 was considered statistically significant.

Statistical analysis

Single fruit weight and soluble solid content were expressed as the mean ± standard error from 10 independent biological replicates by SPSS (V24.0, IBM Corporation, Armonk, NY, USA). These were subjected to one-way analysis of variance (ANOVA), followed by Duncan’s Multiple Range post-hoc test, and the significance level was set to P < 0.01.

Validation of transcripts by quantitative real-time PCR (qRT-PCR)

The expression levels of a set of randomly selected 13 DEGs were validated by a qRT-PCR assay. Total RNA used for RNA-seq was treated with RNase-free DNase I (New England Biolabs, Ipswich, MA, USA) to eradicate all contaminating DNA. A total of 1,000 ng RNA was used for the reverse transcription with PrimeScript™1st stand cDNA Synthesis Kit. qRT-PCR was performed with SYBR Premix Ex Taq (TaKaRa, Dalian, China) on ABI Step One RT-PCR system, according to the manufacturer’s instructions (20 µL reaction mix: 1 µL cDNA, 10 μL 2 ×SYBR real-time PCR premixture, 0.4 µL each 10 µM primer, and 8.2 µL distilled water). Three biological replicates with two technical replicates were performed for each sample. The gene IDs and sequences of 13 primers are listed in Table 1. The PCR program was as follows: 95 °C for 5 min, followed by 40 cycles of 95 °C for 15s, and 60 °C for 30s. Relative expression was normalized to the internal control gene GAPDH gene with 2^−ΔΔCT method (Livak & Schmittgen, 2001). Pearson’s correlation was performed using R software (ver. 3.2.4, R Core Team, 2014) to determine the correlation of gene expression between qRT-PCR and transcriptomic data.

Effect of drought stress on physiological traits and antioxidant activities

Two irrigation treatments were applied to pear trees over the course of three continuous years. Irrigated pear trees were well irrigated, whereas pear trees subjected to deficit irrigation were not irrigated over the same period of time. As shown in Fig. 1, rainfalls during the 2018 season were extremely scarce (Fig. 1A), the consequence of which was a severe decrease in single fruit weight and soluble solids content (Figs. 1B and 1C).

RNA-seq and de novo assembly

Paired-end RNA-Seq was performed on 10 cDNA libraries (5 repeats for drought and irrigation). Each sample was independently aligned, processed for quality control, and then normalized. A total of 400,755,040 clean reads (Table S1) were generated, among which more than 71.7% were mapped to the pear genome GCF_000315295.1_Pbr_v1.0_genomic.fna (Table 2). As indicated by FPKM, the expression values showed high correlations (Spearman correlation coefficient (SCC) = 0.99) among biological replicates, which in turn demonstrated that the sequencing quality was satisfactory for subsequent analyses. Principal component analyses (PCA) revealed that the five replicates of each treatment were located nearest to each other (Fig. 2), thereby demonstrating the reliability of our datasets.

Identification of DEGs between field drought and irrigation treatment

In total, 819 DEGs between drought and irrigation were identified by pairwise sample comparisons (Fig. 3A), among which 427 DEGs were upregulated and 392 DEGs were downregulated in comparison to that of irrigation (Table S2). The expression changes of genes in response to field drought are shown in Fig. 3B. The highly expressed (log2FoldChange <−3.5) drought specific genes (Fig. 3C) included gene38569 encoding Probable WRKY TF40, gene 1490 encoding WRKY TF 18, gene 30473 encoding ferritin-4, gene 7768 and gene 6357 encoding 4-hydroxycoumarin synthase 1, gene 5151 encoding histidine-containing phosphotransfer protein 4, gene 16914 encoding protein NIM1-INTERACTING 1, and gene12366 encoding uncharacterized protein LOC103951864 (Table 3). Genes that were highly expressed in irrigated samples but identified in drought samples included gene 27148 encoding GDL79_ARATH GDSL esterase/lipase At5g33370, gene 5286 encoding uncharacterized protein LOC103944059 isoform X1, gene 1170 encoding putative receptor protein kinase ZmPK1, gene13865 encoding gibberellin-regulated protein 11, gene19880 encoding type I inositol 1,4,5-trisphosphate 5-phosphatase CVP2-like isoform X2, and three genes (gene33465, gene39363, and gene34550) encoding palmitoyl-monogalactosyldiacylglycerol delta-7 desaturase (Fig. 3D, Table 4). The specific expression of 2 DEGs, WRKY TF 18 (gene 1490) and NIM1-INTERACTING 1 (gene 16914), was analyzed by RT-qPCR. Consistent with our RNA-seq results, WRKY TF 18 (gene 1490) was highly expressed in drought treatment at a relatively stable expression level, and the transcription of NIM1-INTERACTING 1 (gene 16914) was consistent with the RNA-seq result only in the irrigation (Fig. 4).

Co-expression analysis of DEGs during field drought treatment

In order to investigate the co-expressed genes during field drought stress, all the genes that were differentially expressed between drought and irrigation were statistically clustered into different groups according to their expression profiles. The largest group (Fig. 5A) included 539 genes that predominantly annotated to RLP12_ARATH and increasingly expressed under field drought conditions. Receptor-like protein 12 participated in the perception of CLV3 and CLV3-like peptides to act as extracellular signals regulating meristems maintenance (149/539). ZIFL1_ARATH Protein ZINC INDUCED FACILITATOR-LIKE 1 (120/539), TMVRN_NICGU TMV resistance protein N (90/539), Y3475_ARATH Probable LRR receptor-like serine/threonine-protein kinase At3g47570 (86/539), and WRK40_ARATH Probable WRKY transcription factor 40 were responsible for the regulation of genes responsive to biotic and abiotic stresses (79/539). The second largest group (Fig. 5B) contained 293 genes whose expression increased under field drought conditions. Genes in this cluster were mainly annotated to BAMS_BETPL Beta-amyrin synthase, which catalyzes the formation of the most popular triterpene among higher plants, HDAC6_HUMAN Histone deacetylase 6, HDAC6_HUMAN Histone deacetylase 6, KAP1_ARATH Adenylyl-sulfate kinase 1, chloroplastic, and RAP24_ARATH Ethylene-responsive transcription factor RAP2-4. The third largest group (Fig. 5C) contained 35 genes whose expression decreased under field drought conditions.

Functional analysis of DEGs between drought and irrigation

Functional analysis was performed to locate enriched GO terms and KEGG pathways involving the DEGs. As shown in Table 5, DEGs were significantly assigned to microtubule (GO:0005874), polymeric cytoskeletal fiber (GO:0099513), and supramolecular complex (GO:0099080) in the cell component (CC) category. In the molecular function (MF) category, DEGs were primarily assigned to microtubules motor activity (GO:0003777), motor activity (GO:0003774), and microtubules binding (GO:0008017). In the biological process (BP) category, DEGs were mainly assigned to microtubules-based movement (GO:0007018) and the movement of cell or subcellular components (GO:0006928). These results demonstrate that DEGs involved in binding, transport, and movement were critical during drought stress.

KEGG pathway enrichment analysis revealed that DEGs were notably enriched in plant monoterpenoid biosynthesis, flavonoid biosynthesis, diterpenoid biosynthesis, cysteine and methionine metabolism, phenylpropanoid biosynthesis, and carotenoid biosynthesis (Fig. 6, Table S3), suggesting specific metabolic events during drought. DEGs were identified using the log2 fold change of the transcript level in field drought compared to the irrigation, and were mapped into the related metabolic pathways (Table 6), thereby revealing a significant impact of field drought on secondary metabolism. Field drought modulated the expression of many DEGs that codify for structural enzymes of the monoterpenoid biosynthesis, flavonoid pathway, and phenylpropanoid biosynthesis (Table 6); the majority of these genes were downregulated under field drought. Four DEGs including the salutaridine reductase-like (SalR) gene family (gene18404, gene39888, gene39889, and gene10010) and nerolidol synthase 1-like (gene237) were involved in monoterpenoid biosynthesis, all of which were downregulated (Table 7) in response to drought stress. Drought modulated the expression of the majority of the structural flavonoid genes (Table 7), most notably three 3,5-dihydroxybiphenyl synthase-like (gene7767, gene7762, and gene 6358), one leucoanthocyanidin reductase-like isoform X1 (gene3879), one BAHD acyltransferase At5g47980-like (gene7261), one salutaridinol 7-O-acetyltransferase-like (gene10701), one vinorine synthase-like (gene34704), and 4-hydroxycoumarin synthase 2 (gene7760). All the aforementioned genes were upregulated by drought. The specific expression of 4 DEGs SalR (gene39889), 3,5-dihydroxybiphenyl synthase (gene7767), BAHD acyltransferase (gene7261), and 4-hydroxycoumarin synthase 2 (gene7760) was analyzed by RT-qPCR, and proved consistent with our RNA-seq results of high expression in drought treatment at a relatively stable expression level (Fig. 4).

Differentially expressed transcription factors under drought stress

Transcription factors (TFs) play key regulatory roles in plant signaling responses, those which activate or inhibit gene expression at the transcriptional level in response to stress. Field-drought treatment led to a number of TFs being differentially expressed (Fig. 7). In total, 4438 differentially expressed TFs were identified, belonging to 30 TF families such as bHLHs (basic helix-loop-helix), NAC (NAM/ATAF/CUC), MYB (v-myb avian myeloblastosis viral oncogene homolog), ERF (ethylene-responsive element binding factor), C2H2s and C3Hs (C2H2 and C3H zinc-finger proteins), WRKYs (WRKY proteins), and bZIPs (basic region-leucine zipper).

Validation of DEG-based gene expression

In order to validate the RNA-Seq gene expression results, qRT-PCR was performed to evaluate the expression levels of the 13 randomly selected DEGs in irrigation vs field-drought conditions (Table 7). As shown in Fig. 4, the expression of the 13 DEGs was largely identical between RNA-Seq and qRT-PCR in spite of certain differences in the absolute fold change. The verified results from the qRT-PCR demonstrated trends similar to the transcriptomic results, which suggests that these DEGs could play significant roles in the regulation of production performance under field-drought conditions.