We collected material from M. sieversii from Daxigou of Huocheng (T3) and Xinyuan population (T7) of Xinjiang at the same time in May 2014 (Table S1). Leaves were immediately frozen in liquid nitrogen and stored at −80°C until needed. Based on the optimal cetyltrimethyl ammonium bromide (CTAB) extraction method, total RNA from equal mixed samples (ten individual plants) was extracted using LiCl purification (Meisel et al., 2005). The RNA purity was evaluated with Nanodrop 2,000 Spectrophotometer (Thermo Fisher Scientific, USA). After RNase-free DNase I treatment (New England BioLabs, USA) to remove residual DNA, cDNA library construction and Illumina HiSeq2500 sequencing were performed by Biomarker Technologies Co., Ltd. (Beijing, China).

EBSeq software (an empirical Bayes hierarchical model for inference in RNA-seq experiments) was applied to perform differential expression genes (DEG) analysis. The process utilizes the well-established Benjamini-Hochberg method, which was corrected for significance of p-values to generate adjusted p-values. During screening, the false discovery rate (FDR) was ≤0.01 and the fold change (FC) was ≥2 (Leng et al., 2013). Combined with the transcriptome sequencing results, MsHsp20 candidates were identified by alignment with the M. domestic genome (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

The MEME Suite version 4.11.3 (http://meme-suite.org/tools/meme) was used to confirm conserved domains of the MsHsp20 protein. The full amino acid sequences of Hsp20 members from M. domestica, Pyrus bretschneideri, and Prunus mume were obtained from NCBI. Gene IDs are shown in Table S3. MEGA 6.0 software was used to construct an unrooted neighbor-joining phylogenetic tree based on 1,000 bootstrap test replicates, pair-wise deletion and a Poisson model.

Based on RNA-Seq data, the full-length open reading frame of MsHsp16.9 was PCR-amplified using primers containing the restriction sites NcoI and BstEII (Table S2). PCR products were firstly cloned into the pEASY-T1 vector and sequenced (TransGen Biotech, China). Following enzyme digestion, the MsHsp16.9 fragment was then inserted into the binary plant vector pCAMBIA3301 under the control of the CaMV35S constitutive promoter (Figure S1). Recombinant vectors were subsequently transformed into Agrobacterium tumefaciens strain LBA4404 using standard heat-shock method.

Arabidopsis (Ecotype Columbia) transformation was carried out by the floral dip method (Clough and Bent, 1998). T4 homozygous transgenic progeny lines were obtained through phosphinothricin-resistant screening and molecular identification. The CTAB method was used to extract genomic DNA from young leaves (Ahmed et al., 2009). Positive pCAMBIA3301 transgenic lines were detected based on the sequences of the 35S promoter of pCAMBIA3301 (Table S2). Primer sets used for the p35S::MsHsp16.9 transgenic lines recognized the forward sequence of the 35S promoter and the reverse sequence of MsHsp16.9 (Table S2).

To analyze heat tolerance in Arabidopsis plants, consistent seedlings were chosen for stress treatment. At 5 days, plantlets in solid 1/2 MS medium were subjected to heat shock at 45°C for 3 h. Eight-leaved plantlets and bolting date seedlings were cultured in a growth chamber at 22 ± 2°C with treatment at 45°C for 16 or 48 h (light intensity = 150 μmol m⁻² s⁻¹). All treatments were under a 16:8 h photoperiod. When all plants showed symptoms of severe wilting, plants were moved back to the normal growth environment.

After heat stress, we estimated the survival rate and root length of 5-day plantlets. The activity of superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT) was determined, as was the malondialdehyde (MDA) content, in eight-leaved plantlets using 752-UV spectrophotometry under heat stress and normal conditions (Wei et al., 2016). After recovering normal growth of florescence, Arabidopsis growth indicators including the rosette diameter, stem length, and the size of the silique were measured and analyzed.

Histochemical staining with 3,3-diaminobenzidine (DAB) or nitro-blue tetrazolium (NBT) was performed to analyze the production of H₂O₂ and O₂- (Shi et al., 2010). O₂- was measured as described previously (Wei et al., 2016), and H₂O₂ levels were measured according to the instructions supplied with the H₂O₂ Assay Kit (KeyGEN BioTECH, China).

PCRs contained 10 μl of SYBR I (SYBR Green qPCR kis, Roche) and reactions was performed using iQ5.0 (Bio-Rad, USA) real-time detection system. The reference genes, MsActin and AtActin, were used as endogenous controls for M. sieversii and Arabidopsis, respectively. iQ5.0 Optical System Software version 2.1 was used for collecting the data. The 2-ΔΔCt method was used to calculate relative gene expression levels from three biological replicates.

All data are presented as the means ± standard error (SE) of at least three replicates. The Student's t-test was used to test the significance of differences between the control plants and transgenic lines. Asterisks (^*, ^**, or ^***) indicate a significant difference between the controls and transgenic plants at P < 0.05, 0.01, or 0.001, respectively.

To investigate the genes response to heat stress, we compared the expression pattern of T7 and T3, and the transcriptome data has deposited in the SRA database of NCBI (T3 accession number: SRR6027256; T7 accession number: SRR6027926). Mixed RNA samples were used to generate a cDNA library, yielding 25,641,037 and 28,821,560 clean reads for T7 and T3, respectively. The guanine–cytosine contents (GC) was 47.64 and 47.55%, and two Q30 base percentage was not less than 85.01% (Table S4). 2,518,193 contigs, 26,013 transcripts and 62912 unigenes were identified on the basis of Trinity assembly (Table S5). We analyzed the expression of the unigenes in T7 and T3. A total of 2728 DEGs were identified via clustering analysis (Figure 1) and divided into 25 groups according to the COG classification (Figure 2). In all COG categories, posttranslational modification, protein turnover and chaperones were located in the third place. In order to confirm the validity of the RNA-Seq data, we selected 18 genes at random that were significantly upregulated and downregulated in T7 and T3 for further analysis by quantitative RT-qPCR (Table S2). The relative expression levels of these genes were similar to those determined from the respective RNA-Seq data. A high correlation (R² > 0.7877) was found between the RT-qPCR and RNA-Seq results (Figure 3, Supplementary Data Sheet 1), confirming the accuracy of the RNA-Seq data.

A total of 12 candidate Hsp20 gene sequences were identified from RNA-Seq DEG data, which were annotated using Swiss-Prot (Apweiler et al., 2004), COG (Tatusov et al., 2000), KOG (Koonin et al., 2004), KEGG (Kanehisa et al., 2004), Pfam (Finn et al., 2014) GO (Ashburner et al., 2000), and non-redundant (Nr) annotation databases with BLAST parameters E ≤ 10⁻⁵ and HMMER parameters E ≤ 10⁻¹⁰ (Table S3). The identified genes showed significant differential expression between T7 and T3, and up-regulation in T3 indicated a common response. Among these genes, the log2FC value was greater than 2 for all except c50641.graph_c0, c55233.graph_c0, and c56990.graph_c0 (Table S3). The C61701.graph_c0 gene, encoding a 16.9 kDa class I Hsp, was expressed highly at all times in both T3 andT7 (Table S3).

The amino acid sequence of MsHsp20 proteins ranged from 136 (>c50697.graph_c0) to 243 (>c34205.graph_c0) residues in length, and all possessed a conserved alpha-crystallin domain (ACD) except for >c30087.graph_c0 and >c33770.graph_c0 based on analysis by MEME (Figure 4A). In addition to >c34205.graph_c0, >c55233.graph_c0 and >c52828.graph_c0, all MsHsp20 members share similar motif, motif1, motif 2 and motif 3 regions, and both nucleotide and amino acid sequences are highly conserved (Figures 4B,C). These motifs are conserved in M. sieversii, and although the functions of these motifs are not yet clear, the presence of similar conserved motifs likely reflects common functions.

To further determine the classification characteristics among MsHsp20 proteins, a phylogenetic tree was constructed with well-supported bootstrap values (1000 replicates), which contained MdHsp20 and PbHsp20 full-length protein sequences. This resulted in the identification of six distinct clusters (classes I–VI; Figure 5). Hsp20 proteins from M. domestica, P. bretschneideri, and P. mume are present in all classes, and M. sieversii sequences are found in all classes. c61701.graph_c0 has high conservation and belongs to class I, together with c56990.graph_c0, c62976.graph_c0, c47100.graph_c0, and c50641.graph_c0. MsHsp20 proteins in classes I and III include c52828.graph_c0, c33370.graph_c0, and c50697.graph_c0, while c55233.graph_c0 and c52423.graph_c0 represent class IV and V, respectively. Other members belong to Class VI.

To further verify that MsHsp20 members are involved in plant responses to environmental stresses, we used RT-qPCR to determine the expression profile of 10 MsHsp20 genes in aseptic seedlings (grown at 25 and 42°C) in response to heat shock (Table S6; Figure 6). The expression levels of the 10 MsHsp20 genes appeared rather low at 25°C, but all were significantly up-regulated in response to heat stress at 42°C (Figure 6). This result indicated that these MsHsp20 genes could respond to heat stress. The primers used are listed in Table S2. To visualize the gene expression patterns, we created a heat map in which each line represents genes with significant differential expression. (Figure S2) Compared with other genes, the expression of >61,701.graph_c0 and >50,697.graph_c0 was relatively lower (Figure S2).

To examine the role of MsHsp16.9 in Arabidopsis, we constructed transgenic plants harboring p35S::MsHsp16.9 (Figure S1). We selected pCAMBIA3301 vector control transgenic lines by spraying Basta and molecular identification (Figure 7A). Similarly, 12 positive MsHsp16.9 overexpression lines were initially screened with Basta, and PCR was then performed with specific primers, resulting in the amplification of a 471 bp MsHsp16.9 fragment (Figure 7B). We performed RT-qPCR analysis using cDNA from p35S::MsHsp16.9 lines 3, 6, 8, and 11. Compared with untransformed control plants (WT), expression of MsHsp16.9 was 400–1500-fold higher. We chose relatively high expressing lines (3, 8, and 11) for subsequent experiments (Figure 7C). Tissue localization of MsHsp16.9 expression was then investigated, and expression was particularly high in the leaves and pods, and lower in stems and flowers (Figure 7D).

We observed the growth of MsHsp16.9 transgenic Arabidopsis plants and compared them with WT plants and pCAMBIA3301 vector control plants in 1/2 MS medium. Three overexpression lines showed no obvious changes compared with the control (Figures 8A, 9B). To investigate whether MsHsp16.9 expression was associated with heat tolerance, we subjected transgenic plantlets together with WT and pCAMBIA3301 transgenic plants possessing 2–4 leaves to heat at 42°C for 3 h (Figures 8A, 9B). The survival rate of the three lines of transgenic plants was 16, 25, and 41%. However, survival of WT and pCAMBIA3301 transgenic plants was only 8% (Figure 8C). In addition, the average root length of transgenic plants was longer than control plants (Figure 8D). These phenotypic differences suggest overexpression of MsHsp16.9 may weaken the plant sensitivity to heat stress.

To investigate whether heterogenous expression of the MsHSP16.9 gene alters transgenic plants heat responses, transgenic and control Arabidopsis with 8 to 10 leaves were incubated at 45°C for 16 h (Figure 9A). The results showed that transgenic lines were less sensitivity than WT and VC control plants, and displayed less withering, indicating more stable turgor pressure and slower onset of senescence (Figure 9B). After 7 days recovery from the heat treatment, transgenic plants exhibit varying degrees of growth and some leaf death. By contrast, almost all WT and VC plants that developed leaves had died (Figure 9C). At 10 days after heat treatment, surviving WT and VC plants had begun to grow again, whereas most transgenic lines had completely restored growth by this stage and were close to the bolting stage (Figure 9D). The plant survival rate was further analyzed, and the results showed that both control and transgenic plants under normal growth conditions had a similar growth rate, but control plants was less than that of transgenic plants under heat stress (Figure 9E). Additionally, we investigated oxidase system activity and membrane damage. The activity of SOD, POD and CAT was clearly increased in transgenic plants (Figures 9G–I). Consistent with this, elevated MDA levels were found in control plants under stress conditions, compared with the lower levels observed in transgenic lines (Figure 9F).

To further assess the effects of heat treatment on reproductive growth, WT, VC, and MsHsp16.9-overexpressing transgenic bolting stage plants were incubated at 45°C for 48 h then recovered until the full blossom stage (Figures 10A–C). Rosette leaf diameter, stem length, and pod length were measured in WT, VC and transgenic seedlings after recovery. Regarding the overall plant structure, control and transgenic plants showed different degrees of bleaching and flowering time delay. Furthermore, transgenic plants had formed pods, while control plants had lost nearly all reproductive ability (Figures 10F,K). The main reason is that the stem and branches of control plants adopted a dwarf phenotype, leading to the failure of flowering and pollination (Figures 10D,E,G,J). Consistently, the smaller rosette leaves and single blades of WT and VC plants also had an adverse effect on blossoming (Figures 10H,I).

Heat stress can increase the generation of reactive oxygen species (ROS), thus we investigated whether transgenic plants accumulated less ROS under heat stress treatment by examining O₂- and H₂O₂ accumulation using NBT and DAB staining, respectively. Under normal growth conditions, DAB and NBT staining was minimal in transgenic and non-transgenic plants (Figures 11A,C). However, staining was increased following heat treatment in all plants, but injury in control plants was clearly more serious than in transgenic plants (Figures 11B,D). Regarding O₂- and H₂O₂, levels of both were higher in WT and VC plants (Figures 11E,F). These results further confirmed the protective role of MsHSP16.9 in heat stress tolerance.

The plant hormone abscisic acid (ABA) is a major regulator of plant growth and a variety of stress environments (Finkelstein et al., 2002; Xiong et al., 2002; Himmelbach et al., 2003). To assess whether the role of MsHsp16.9 in abiotic stress responses involves ABA, we compared sensitivity to exogenous ABA during germinative and post-germinative growth in WT and transgenic plants. The results showed that MsHsp16.9-overexpressing plants did not exhibit an obvious difference in their response to exogenous ABA during germination or growth (Figure 12A). However, the expression of genes in ABA-mediated signaling pathways was distinctly changed compared with WT plants. The 9-cis epoxycarotenoid dioxygenase (NCED) encoded by the AtNCED3 gene can catalyze the limiting step of ABA biosynthesis (Roychoudhury et al., 2013). AtNCED3, together with ABA-responsive gene, AtABI, AtABI5, AtSNF4, and AtABER2, was up-regulated in all transgenic lines compared with the controls under normal and ABA treat conditions (Figure 12B), especially in 35S:MsHsp16.9-L11 (TP11). This difference in ABA biosynthetic pathway gene expression between MsHsp16.9 transgenic and control plants suggests that heterogenous expression of MsHsp16.9 may increase ABA biosynthesis and accumulation in plants, and raises the possibility that MsHsp16.9 acts as a positive regulator in ABA signaling in Arabidopsis.

To explain the phenotypic characteristic of the MsHsp16.9 transgenic plants under heat stress and have an insight into the molecular mechanisms of MsHsp16.9-mediated regulation of the heat stress response, we focused on heat stress-related genes that are differentially expressed between WT and MsHsp16.9 transgenic plants. These genes were divided into two groups: regulatory factors and functional proteins. The regulatory proteins mainly included heat-shock factor (HSF) family members (HSF1D, HSF1E, HSFA3, and HSFA4A) as well as the dehydration-responsive element binding protein (DREB2A). All HSF genes were up-regulated in TP11 but not in other lines. The expression of Hsp70 was also consistent with HSF. These results, together with the observed phenotypes and high expression of MsHsp16.9 in TP11, suggests that MsHsp16.9 may activate Hsp70 and HSF to maintain protein stability, directly or indirectly protecting cells from damage during high expression. Besides, DREB2A has pivotal effect on regulating the heat response in all transgenic lines including 35S:MsHsp16.9-L3 (TP3), 35S:MsHsp16.9-L8 (TP8) and TP11 (Figure 13). Regarding functional proteins, the high expression of CAT and ascorbate peroxidase (APX) genes in transgenic lines TP3 and TP11 may explain the observed alleviation of oxidative damage (Figure 13). In addition, arginine decarboxylase (ADC1), S-adenosylmethiomine decarboxylase (SAMDC) and 1-pyrroline-5-carboxylate synthase (P5CS), all associated with secondary metabolism, were up-regulated in transgenic lines under heat stress compared with WT and VC plants (Figure 13). Similar results were observed for the early response to delydration Stress (ERD) gene. These osmotic regulation-associated proteins was crucial for cellular homeostasis and be affected by MsHsp16.9 under heat stress.

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.