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

Drought acts as a major abiotic stress that hinders plant growth and crop productivity. It is critical, as such, to discern the molecular response of plants to drought in order to enhance agricultural yields under droughts as they occur with increasing frequency. Pear trees are among the most crucial deciduous fruit trees worldwide, and yet the molecular mechanisms of drought tolerance in field-grown pear remain unclear. In this study, we analyzed the differences in transcriptome profiles of pear leaves, branches, and young fruits in irrigation vs field-drought conditions over the growing seasons. In total, 819 differentially expressed genes (DEGs) controlling drought response were identified, among which 427 DEGs were upregulated and 392 DEGs were downregulated. Drought responsive genes were enriched significantly in monoterpenoid biosynthesis, flavonoid biosynthesis, and diterpenoid biosynthesis. Fourteen phenylpropanoid, five flavonoid, and four monoterpenoid structural genes were modulated by field drought stress, thereby indicating the transcriptional regulation of these metabolic pathways in fruit exposed to drought. A total of 4,438 transcription factors (TFs) belonging to 30 TF families were differentially expressed between drought and irrigation, and such findings signal valuable information on transcriptome changes in response to drought. Our study revealed that pear trees react to drought by modulating several secondary metabolic pathways, particularly by stimulating the production of phenylpropanoids as well as volatile organic compounds like monoterpenes. Our findings are of practical importance for agricultural breeding programs, while the resulting data is a resource for improving drought tolerance through genetic engineering of non-model, but economically important, perennial plants.


INTRODUCTION
Along with an increasing global population, drought is becoming one of the most persistent factors that limits agricultural production and food security around the world, especially in arid and semi-arid regions (Mittler, 2006). In turn, drought is responsible for losses in the multibillions of dollars annually (Fahad et al., 2017;Lesk, Rowhani & Ramankutty, 2016). China is facing a perilous water crisis in which 50% of the national territory is located in arid and semi-arid regions (Hu & Zhang, 2001). The temporal-spatial distribution of annual precipitation causes 26.7% of the national land territorial area in Northwest China to have arid and semiarid climates, a region where drought is common. Predicting drought severity is difficult, and to do so requires consideration of several factors such as rainfall amount and distribution, evaporative demands, and the moisture storing ability of soils (Saud et al., 2017;Tadesse & Melkam, 2016). Globally, several management strategies have been implemented for improved crop production under drought environments (Bodner, Nakhforoosh & Kaul, 2015;Fahad et al., 2017). Among these, the development of crop varieties with an increased tolerance to drought functions as an important and effective strategy to combat drought.
Plants cope with water deficiency by complex mechanisms from molecular, biochemical and physiological processes at the cellular or whole plant level (Bray, 1997;Goufo et al., 2017;Huber & Bauerle, 2016;Ruggiero et al., 2017;Zandalinas et al., 2017). With the advent of new high throughput ''-omics'' technologies like proteomics and transcriptomics, notable strides have been made towards understanding the molecular mechanisms that regulate tolerance to drought. Previous studies have demonstrated signal transduction of drought stress perception to the nucleus via complex cellular signaling networks involving second messengers. These include reactive oxygen intermediates (ROIs) and calcium, calciumassociated proteins, and kinase cascade such as mitogen-activated protein (MAP) (Bray, 1997;Chen et al., 2002;Huber & Bauerle, 2016;Knight & Knight, 2001;Liu et al., 1998;Zandalinas et al., 2017). Drought stress signaling cascades are comprised of many stressresponsive genes. These include molecular chaperones such as late embryogenesis abundant (LEA) proteins and heat shock proteins (HSPs) that function as effector molecules. Other examples include transcription factors (TFs) like members of the APETALA2/ethyleneresponsive element binding protein (AP2/EREPB), a basic leucine zipper (bZIP), WRKY, and MYB proteins that act as regulator molecules (Shinozaki & Yamaguchi-Shinozaki, 2007;Song et al., 2005;Wang, Vinocur & Altman, 2003). The physiological and molecular mechanisms of plant responses to drought have been extensively studied in model plants with dehydration treatments in controlled laboratory or greenhouse conditions (Li, Xu & Huang, 2016;Wang et al., 2018;Zarafshar et al., 2014). However, results from these studies most often translate poorly to field-grown plants. Clarifying the molecular mechanisms that regulate drought tolerance from crops grown under field conditions will facilitate a more thorough grasp of the complex interactions between drought response and environmental factors that crops encounter in the field during the growing season. As such, the task of developing an improved understanding of molecular elements in responsiveness to field drought in non-model plants will aid in both traditional and modern breeding applications towards improving stress tolerance.
Pear is one of the most vital fruit crops in the world and the second major crop among deciduous fruits in China after apples (Silva et al., 2014). The crop has considerable value both economically and in terms of personal health. In China, pear is primarily grown in the Northwestern region, accounting for 60 percent of pear production in the country. YuluXiangli (Pyrus spp) is an improved pear cultivar that is highly tolerant to drought, and it is an ideal source for examining genomic responses to drought in order to explore valuable tolerance genes (Okubo & Sakuratani, 2000;Zong et al., 2014). The full genome sequencing and resequencing of multiple pear cultivars (Huang et al., 2015;Li, Xu & Huang, 2016;Wang et al., 2018;Wu et al., 2013) have enabled several transcriptome studies of drought responses in pear, thereby revealing a broad, multifaceted response to drought. Such a response features coordination between phytohormone signaling pathways, the reduction of photosynthetic gene expression, and the alteration in expression of genes involved in stress-induced leaf senescence. These studies, however, have been restricted to greenhouses under certain durations of drought stimuli treatment as opposed to field conditions that use early time points with samples exclusively from leaves (Li, Xu & Huang, 2016;Wang et al., 2018).
The primary objectives of the present study were to identify differentially expressed genes (DEGs) and to compare the gene expression patterns in leaf, branch, and fruit tissue of pear in response to drought induced by withdrawal of irrigation in the field. The findings will provide an unrivaled resource for understanding the mechanisms underlying drought resistance in pear.

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 log2transformed 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 TM 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 Table 1 The gene IDs and primer sequences for qRT-PCR.

ID
Primer 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.

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   (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 Ethyleneresponsive 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 4hydroxycoumarin 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 fielddrought 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.   drought tolerance through programmed selection with precise strategies of stress-testing, particularly in light of ongoing global climate change. In the present study, we identified differentially expressed genes under field-drought stress and irrigation control with RNA-Seq in the pear cultivar YuluXiangli. A total of 819 DEGs were detected, and 4,438 TFs were differentially expressed between drought and irrigation control. Our findings represent valuable information on transcriptome changes in response to drought. Drought responsive genes are mainly enriched in biosynthesis-related pathways-monoterpenoid biosynthesis, flavonoid biosynthesis, and diterpenoid biosynthesis-and they belong mainly to bHLHs, NAC, MYB, ERF, C2H2s, and C3Hs, as well as to WRKYs transcription factor families. Our analysis provides a solid foundation for both the identification and the functional analysis of potential candidate genes related to drought tolerance. The prolonged and severe field drought imposed in this experiment modulated the accumulation of phenylpropanoids, flavonoids, monoterpenoid biosynthesis, and several volatile organic compounds in the pear. Previous studies demonstrated the droughtmodulated accumulation of phenylpropanoids, flavonoids, terpenoids, and carotenoids under drought (Li et al., 2018;Murphy & Zerbe, 2020;Savoi et al., 2016;Sircelj et al., 2005). This accumulation acted as antioxidants and protected plants from the adverse effects of drought conditions (Nichols, Hofmann & Williams, 2015). Our study demonstrated modulation of the biosynthetic pathways of phenylpropanoids and flavonoids by drought stress at the transcript level, leading to an enhanced accumulation of derivatives of benzoic and cinnamic acids as well as several flavonoids. This was congruent with previous results (Li et al., 2018;Murphy & Zerbe, 2020;Nichols, Hofmann & Williams, 2015;Savoi et al., 2016;Sircelj et al., 2005). Five of 14 phenylpropanoid DEGs as well as all of the flavonoid DEGs were upregulated under drought stress, the result of which enhanced the concentration of accumulation within these compounds. Flavonoid aggregation in cytoplasm is capable of effectively detoxifying drought-induced harmful H 2 O 2 molecules. In the present study, the elevated flavonoid aggregation was induced by drought stress condition, supporting previous results in Achillea pachycephala Rech.f. (Gharibi et al., 2019), Brassica napus (Rezayian, Niknam & Ebrahimzadeh, 2018), Arabidopsis (Nakabayashi et al., 2014), grape (Degu et al., 2015Savoi et al., 2016), and white clover (Ballizany et al., 2012). The physiological and molecular mechanisms underlying the drought-induced accumulation of these compounds to modulate phenylpropanoid as well as the flavonoid biosynthetic pathway need to be further elucidated by integrated transcriptome and metabolite profiling.
The monoterpenoid biosynthesis was significantly modulated by the prolonged and severe field drought conditions in the present experiment. Plant terpenes were synthesized in the plastids through the 2C-methyl-D-erythritol-4-phosphate pathway (MEP), and in the cytosol through the mevalonate (MVA) (Tholl, 2006). A number of terpenoid metabolites were involved in adaptation to adverse environments (Pichersky & Raguso, 2018;Murphy & Zerbe, 2020), including biotic and abiotic stresses; however, the knowledge of drought-modulated regulatory mechanism of monoterpene biosynthesis is limited (Zhang et al., 2019). All of the four-terpene synthase (TPS) genes encoding salutaridine reductase (SalR) and nerolidol synthase 1 involved in monoterpene biosynthetic pathway were downregulated under drought conditions. Salutaridine reductase catalyzes the stereo specific reduction of salutaridine to 7(S)-salutaridinol, nerolidol synthase 1 converts geranyl diphosphate (GPP) into S-linalool, and farnesyl diphosphate (FPP) into (3S)-E-nerolidol in the biosynthesis of morphin (Ziegler et al., 2009). Morphine resides within the diverse class of metabolites called benzylisoquinoline alkaloid, and drought stress, it has been noted, can increase alkaloids in opium poppy (papaver soniniferum) (Szabó et al., 2003). Our results were unlike previous findings in several plants (Selmar & Kleinwächter, 2013), those such as Chrysopogon zizanioides (Ziegler et al., 2009) and grapevine (Griesser et al., 2015;Savoi et al., 2016). Ziegler et al. (2009) reported upregulation of the gene encoding Salutaridine reductase under drought stress specifically in leaf tissue of Chrysopogon zizanioides. Drought-induced monoterpene production was observed in several plants (Selmar & Kleinwächter, 2013) including grapevine leaves (Griesser et al., 2015;Savoi et al., 2016). Six TPS genes, one of which included the nerolidol synthase 1-like gene, were differentially expressed in response to abiotic stresses in Santalum album (Zhang et al., 2019). Further biochemical and transcriptomic profiling is needed to address terpenoid biosynthetic pathways and their spatiotemporal regulation in response to adverse drought stress.
Transcription factors (TFs) modulate diverse transcriptional regulation and play significant regulatory roles in plant signaling responses to developmental and environmental changes (Ulker & Somssich, 2004;Osakabe et al., 2014;Yu et al., 2015). In the present study, 4,438 differentially expressed TFs were identified to promote or suppress abiotic stress responses, including the bHLHs, NAC, MYB, ERF, C2H2s, C3Hs, and WRKY families. WRKY TFs have been reported to be involved in drought stress responses through the ABA signaling pathway (Ulker & Somssich, 2004;Osakabe et al., 2014). Overexpression of ZmWRKY58 enhances the drought and salt tolerance in transgenic rice (Cai et al., 2014). Drought-responsive WRKY TFs TaWRKY33 and TaWRKY1 confer the transgenic Arabidopsis plants drought and/or heat resistance (He et al., 2016). The cotton WRKY TF GhWRKY33 reduces transgenic Arabidopsis resistance to drought stress (Wang et al., 2019). In the present study, a total of 79 WRKY genes induced by field drought treatment were grouped in cluster 1, the majority of which were upregulated. Li, Xu & Huang (2016) reported 637 transcription factors responsive to dehydration in pear, among which 45 WRKY genes were differentially expressed. Huang et al. (2015) classified a total of 103 WRKY TFs in the pear genome, demonstrating an improvement of tolerance to drought by manipulating the PbWRKYs. Therefore, WRKY TFs may play significant roles in regulating drought stress responses.

CONCLUSION
We utilized deep sequencing technology to investigate the transcriptome profiles in pear leaves, branches, and young fruits in response to the prolonged field drought induced by irrigation withdrawal. A total of 819 DEGs were detected, and 4,438 TFs were differentially expressed between drought and irrigation control, presenting valuable information on transcriptome changes in response to drought. We illustrated the flavonoids and monoterpenoid biosynthesis-related genes specifically expressed in drought and irrigation control during field-grown season in pear. Validation of gene expression by 13 randomly selected genes was in correspondence with transcriptomic results. Several candidate genes including flavonoid and terpenoid genes, transcription factors, and drought-responsive elements, were involved in transcriptional regulation of plant response to drought. Such information is important to germplasm management and in endeavoring to improve pear productivity.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This work was supported by the Shanxi Province Natural Science Foundation (201801D121255, 20210302124122), the Research Subject of Agricultural Science and Technology Innovation of Shanxi Academy of Agricultural Sciences (YCX2018D2YS14), the Project of Scientific and Technological Innovation research of Shanxi Academy of Agricultural Sciences (YCX2020SJ10), and the China Agriculture Research System (CARS-28-28). There was no additional external funding received for this study. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.