ABBREVIATIONS

INTRODUCTION

Heat stress is one of the most severe challenges for the world food security [1]. It was shown that the average global temperature is increasingly rising in the past decades [2]. The increased emission of greenhouse gas, e.g., carbon dioxide, methane, nitrous oxides, production is considered as the major reasons for the elevated temperature [3]. High temperature can reduce growth rate and biomass accumulation in plants [4,5]. This negative effect can be attributed to alterations at the morphological, anatomical, and physiological levels. For example, in sugar cane (Saccharum officinarum), shortened internodes, early senescence and more but smaller tillers were observed to be associated with heat [6]. In addition, reduced cell size, closure of stomata and presumably lowered photosynthesis were also found to occur in various plant species in response to heat [7,8]. Heat stress can harm plants at both vegetative and reproductive stages. Under extreme circumstances, it can lead to plant death due to the irreversible damage to cellular components, increased respiration and synthesis of reactive oxygen species, disruption of cell membranes, destruction of structure and functions of proteins and interference of important metabolic pathways [9,10].

Plants are sessile organisms. Therefore, they need to develop various physiological and molecular mechanisms to assist them for the adaptation to elevated atmospheric temperature. Plants might adopt avoidance mechanism to cope with heat stress which includes improvement of cooling via elevated transpiration, alteration of leaf angle, optimization of cell membrane lipid composite, and stronger rooting [11]. Even plants develop mechanisms to help them avoid the negative effect under heat, high temperature nevertheless can exceed the ability of the avoidance ability in plants under many circumstances. Another mechanism that plants develop is to adjust their physiological metabolisms so they can tolerate the adverse effect. This includes increased content of compatible compounds, activated free radical scavengers, increased abundance of chaperon proteins, ion transporters, late embryogenesis abundant (LEA) proteins and accelerated transcription of transcription factors encoding heat shock proteins [3]. A third mechanism that some plants develop is the escape strategy with which some plants fulfill their life cycle faster with a comprised yield [12].

The elevated temperature can induce alterations in physiological, morphological, and biochemical levels via changing expression levels of multiple genes and pathways. The recent advancement of genomic tools such as RNA-Seq greatly accelerates the understanding of the molecular responses at the whole genome level. In addition, the development of a few informatics database such as GO, COG, and KEGG allow the identification of the most enriched pathways in response to stress conditions. A number of transcriptomic or proteomic profiling studies have shown that hundreds of genes or proteins are involved in response to high temperature in soybean. However, these studies were focused on either the responses in roots or vegetative stage of the soybean plants [13–15]. Although high temperature can cause growth reduction and yield loss for plants at all development stages, it has been found that reproductive stage plants are more sensitive to the stress. Our aim of this study was to use a high-throughput RNA-Seq profiling technique to study the molecular mechanisms in reproductive stage soybean in response to heat. The compromised physiological processes and the diverse genes and pathways identified in this study indicated that heat causes severe adverse effects to soybean and the stressed plants employ complicated mechanisms to cope with this stress. The dissection of this process will not only elucidate the basic physiological mechanisms associated with this stress, but it also provides targets for either genetic engineering or molecular breeding for crop improvement.

MATERIALS AND METHODS

Williams-82 soybean (Glycine max L. cv. Williams) was used in this study. The rhizobia were inoculated on the soybean seeds to improve nodule formation following the instruction provided by the manufacturer (Advanced Biological Marketing, OH, USA). In brief, approximately 500 soybean seeds were first dampened with tap water and kept in a container. The seeds were then mixed with the rhizobia inoculant at the rate of 3.0 g/kg of seeds. When the surface of soybean seeds coated with the bacteria inoculant was dried, three seeds were planted in pots filled with propagation mix (Sun Gro^® Horticulture, Agawam, MA, USA) and in total twenty pots were prepared.

When the soybean seeds were germinated, only one seedling was kept in each pot. These plants were continuously grown in a growth chamber with the day/night temperature at 28/25 °C and light intensity at approximately 400–500 µmol·m⁻²·s⁻¹ on a 14 h:10 h light regime. Heat stress treatment was initiated at approximately 45th days (R2 stage) post planting by moving half of the plants into a second chamber with similar settings as the one for their initial growth except a higher temperature setting at 38/25 °C.

Three days after the plants were treated with the high temperature, a Li-Cor 6400 (LI-COR, Lincoln, NE, USA) photosynthesis apparatus was used to measure net photosynthetic rate, transpiration, and stomatal conductance on the top three fully expanded leaves according to the instruction from the manufacturer. The temperatures on detached leaves were measured with a “Flir-One” system (FLIR, Nashua, NH, USA). Relative water content was calculated as:

After the leaves were sampled and weighed (fresh weight), they were immersed for 24 h in distilled water, and weighed again (hydrated weight). The leaves were then dried for 24 h in a 75 °C oven and weighed (dry weight) [16].

The top three fully expanded leaves used for the photosynthesis measurements were frozen in liquid nitrogen immediately after the measurements were finished. These leaves were stored in −80 °C freezer until the analysis of antioxidant enzymes, cDNA synthesis and preparation of RNA-Seq libraries were conducted. The activities of three antioxidant enzymes, i.e., POD, CAT, and SOD, were measured. The detailed procedures were described in our recent publication [17]. Three replicates were used for either control or heat stressed soybean plants. The total protein content was determined using Bradford’s method [18].

We collected three top leaves from each of the control or heat stressed soybean plants and then pooled them as one sample, with 3 replicates. Each sample contains 100 g of leaf tissue and grounded in −80 °C liquid nitrogen. Total RNA was extracted by E.Z.N.A plant RNA kit (Omega Bio-Tek Company, Norcross, GA, USA). Extracted total RNA were measured by NanoDrop spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE, USA) to check the quality and quantity of RNA. Reverse transcriptions were conducted on these RNAs and six libraries were constructed from these cDNAs following the instruction from the manufacture (Illumina Company, San Diego, CA, USA). A QC was conducted on these samples before RNA Seq was performed.

The cDNA libraries were sequenced on a HiSeq4000 (Illumina) sequencing platform using 51 bases, single-end sequencing, and the quality of RNA-Seq data was determined using FASTQC (v0.10.1) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), which was followed by mapping the reads to reference genome by TopHat2 [19], with up to two mismatches. Differentially expressed genes (DEGs) were identified using Cuffdiff following normalization of transcript count information to RPKM (reads per kilobase of exon model per million mapped reads) [20] and genes with a P-value < 0.05 and fold change 2 cutoff were considered to be differentially expressed.

To compare the activities for antioxidant enzymes and the rates for photosynthesis, relative water content, stomatal conductance, three replicates were used for each of these measurements and a student-t test and a 5% P value were used for significance calculation.

RESULTS AND DISCUSSION

The initial molecular alterations under stress conditions in a plant will ultimately lead to re-programming of multiple physiological processes for adaptations. Before dissecting the molecular responses associated with heat in soybean, it is important to assure that some of the typical heat stress associated physiological alterations have been initiated. A number of physiological parameters were measured on soybean plants (R2 stage) that were grown either at control or heat stressed condition. It is not surprising to notice that the temperature of leaves on well-watered soybean plants grown at 38 °C was higher than the temperature of leaves on well-watered soybean plants grown at normal temperature (28 °C) (Figure 1). However, when being compared with our recent study in which soybean plants were grown in drought condition, the increase of temperature under heat is much smaller than that under drought [17]. This indicated that the soybean plants grown under heat condition were actively undergoing transpiration to cool the leaves.

Compared to the well-watered plants, the photosynthesis and relative water content in the heat stressed plants was reduced by 16.7% and 5.3%, respectively (Figure 2). Previous studies have shown that under high temperature certain photosynthesis related proteins can be denatured [21,22]. The slightly increased stomata conductance and transpiration indicated that these soybean plants were experiencing heat stress and evaporating water to lower the leaf temperature (Figure 2).

FIGURE 1

Figure 1. Leaf temperature of soybean plants subjected to CK (control) and H (heat) stress. Leaf has slightly higher temperature under heat stress (38 °C) than control condition (28 °C). The leaf temperature was measured using a “Flir-One” thermal imaging system (Flir, Nashua, NH, USA).

FIGURE 2

Figure 2. Measurement of physiological parameters in plants subjected to control (CK) and heat (H) stress conditions. “a” and “b” indicate significance using a two-sample t test assuming unequal variance. Bar represents SE (N = 3). Significance was defined as P ≤ 0.05.

Reactive oxygen species (ROS) are signaling molecules in plants. However, high concentration of ROS has negative effect to the metabolism and cell structure in plants [23]. Therefore, plants have evolved mechanisms to accurately control equilibrium of ROS in plant cells [24]. We measured the activities of three antioxidant enzymes. Different from our discoveries in drought stressed soybean plants, in which the enzymatic activities for both CAT and SOD were increased [17], in the heat stressed soybean plant, the activities of these two enzymes were not changed. However, the activity of POD was increased (Figure 3). The increase of POD activities in response to high temperature had been discovered in multiple plant species [25,26]. The findings from our study indicated that, under our treatment condition for the reproductive stage of soybean plants, drought and heat stress use different antioxidant enzymes to cope with the over accumulated ROS.

FIGURE 3

Figure 3. Antioxidant enzyme activities in plants subjected to control (CK) and heat (H) stress conditions. A student-T test and a 5% P value were used for significance calculation. “a” and “b” above the column for CK and H indicated the measurements were statistically different. Bar represents SE. For each of the measurements, three replicates were used.

A total of 334.7 Mbp reads (raw reads) derived from leaves collected from both control and heat stressed soybean were generated from the Illumina Hi-Seq 4000 platform. After processing and filtering of raw reads, a total of 161.3 Mbp (clean reads) were obtained. Over 50% of reads from each sample could be mapped to the assembled transcriptome (Table 1).

TABLE 1

Table 1. Classification of total transcript produced in soybean leaves under control and heat stress.

Single reads from the sequencing were analyzed and all the clean reads were compared to the genome sequence by the TopHat2 system. For functional annotation analysis, we used COG, GO and KEGG database with BLASTX method, and obtained 247, 442 and 227 respectively, unigenes that were annotated to these databases (Table 2).

TABLE 2

Table 2. Heat-stress-responsive DEGs annotated to different databases.

FIGURE 4

Figure 4. The degree of expression change for heat stress responsive DEGs shown in a volcano plot. X-axis represents fold change; Y-axis represents significance. Each dot represents a DEG.

A total of 623 unigenes were annotated to all the databases combined. To increase the probability of authentic gene identification, only genes with at least 2-fold changes (FC) were selected for further investigation. A total of 633 annotated genes were found to be differentially expressed in heat stressed plants in which the expression of 417 genes were up-regulated and 216 were down-regulated (FC > 2, P-value < 0.05). The degree of expression change of these DEGs is shown in volcano plot (Figure 4). In our previous drought stress transcriptomic analysis on soybean plants at similar development stage, a total of 2771 annotated genes were found to be differentially expressed [24], indicating the less severe effect on molecular alterations to heat than that to drought.

To get an overview of the functional classification of the assembled unigenes derived from the heat-stressed soybean leaves, GO analysis was performed. A total of 422 DEGs were assigned by KOBAS2.0 and HMMER into 3 broad categories: the biological process, cellular component and molecular function. Under heat stress, “cellular process”, “metabolic process” and “single-organism process” are the highest responsive groups in the biological process category; “cell part” and “cell” groups are the highest groups in the cellular component category; “binding” and “catalytic activity” are the highest groups in the molecular function category (Figure 5).

FIGURE 5

Figure 5. Ontology classification of DEGs in the Heat-stressed soybean. The X-axis represents the functional categories; Y-axis (left) represents percentage of genes falling into this category; Y-axis (right) represents the absolute number of genes in the corresponding categories.

To understand the functionalities of these genes from the perspective of phylogenetic classification of proteins, COG analysis was performed and 247 genes were assigned into 26 functional categories. The most enriched groups were “general function prediction” followed by “signal transduction”, “transcription”, “replication, recombination, and repair” and “carbohydrate transport and metabolism” (Figure 6).

FIGURE 6

Figure 6. Distribution of the DEGs in response to heat stress. X-axis represents the gene functional categories; Y-axis represents the absolute number of genes in the corresponding category.

Similar to what we conducted on genes identified from drought stressed soybean plants [17], a KEGG pathway enrichment analysis was initially conducted to identify enriched pathways in soybean grown under heat stressed condition. To our surprise, only 142 out of 633 DEGs could be recognized by the software and they were assigned to 78 pathways (Supplementary Table S1, Supplementary Figure S1). The limited number of genes on most of these pathways hindered an insightful elucidation of the altered physiological pathways. Therefore, to get a more comprehensive understanding on the major categories of genes differentially expressed under heat, we went through literature for all the differentially expressed genes and grouped them based on their involvement in plant physiological processes. Here we chose a few representative groups to illustrate the mechanisms that soybean plants use to respond to heat condition.

The flowering processes can be either induced or repressed in plants grown under heat condition. In Arabidopsis (Arabidopsis thaliana), elevated temperature stimulates flowering [27], while in chrysanthemum, higher temperature delays flowering [28]. In some cereal plants, elevated temperature delay flower development under short days but promote flowering under long days [29]. In this study, we found that the expression levels of ELF3, GI, FKF1, PRR5, and TOC1 were increased; while the expressions of LHY and FT were decreased (Table 3).

TABLE 3

Table 3. Differentially expressed genes related to flowering regulation.

In addition, the expression of CO differed among different members in this gene family. One common characteristic about all of these genes is that they are controlled by circadian rhythms (Figure 7). Previous studies have shown that ELF3, GI, FKF1, PRR5, TOC1, CO, and FT are positive regulators for flowering, while LHY is a negative regulator for flowering [30–34]. For example, overexpression of CO or FT driven by either a constitutive promoter (35S) or a phloem specific promoter (SUC2) in plants had led to early flowering [35–37]. In contrary, T-DNA lines with abolished FT transcript or transgenic plants using cell ablation method to reduce both CO and FT had delayed flowering [38]. As another example, early flowering was found in LHY-mutant Arabidopsis plants even under non-inductive conditions [39]. While the expression tendency for the majority of these genes indicated the soybean tended to accelerate flowering process under heat condition, it is puzzling to see the lowered expression levels of FT and a member of the CO transcripts as the reduced expression of these genes should delay flowering. An observation of this kind indicated the limitation of using sole transcriptomic analysis to interpret the physiological alterations in soybean in response to heat since post-transcriptional and translational regulations may also be involved in response to heat. However, these findings did demonstrate that the adjustment of flowering time in soybean is one of the mechanisms to cope with heat stress. A combination of transcriptomic and proteomic analysis in the future will provide more unambiguous understanding to the involvement of flowering time regulation in soybean under heat stress.

FIGURE 7

Figure 7. The pathway of “Circadian Rhythm” under heat stress condition. The downregulated genes involved in “Circadian Rhythm” are marked in green; the upregulation genes are marked in red. Reproduced with permission from Kanehisa Laboratories.

Reactive oxygen species (ROS) are tended to be generated in response to stress conditions such as drought and heat [24]. An inhibition of the synthesis or a timely removal of the over-produced ROS compounds from plant cells are needed to prevent delirious reactions with certain cellular components and harmful effect to essential cell structures such as plasma membranes. In this study, we found that among the activities of the three measured antioxidant enzymes, i.e., SOD, POD, CAT, only that for POD was increased (Figure 3). This is different from the responses under drought condition where the activities of SOD and CAT were increased [17]. The differential responses indicated other antioxidant mechanisms may exist for the removal of ROS in heat stressed soybean plants. In this study, the expression levels of a total of 64 of the pentatricopeptide repeat-containing protein (PPR) gene family were found to be increased (Table 4). Members of the PPR gene family is characterized by the harboring of 9 to 15 arrays of 35 amino acid long pentatricopeptide repeats [40]. Multiple physiological functions have been associated with this gene family. For example, it is known that PPR proteins regulate RNA metabolism, such as RNA splicing, editing, and translation [41]. In addition, recent research evidence showed that genes in this family are also involved in electron transport, reactive oxygen species generation and abiotic stress resistance. Different from the protective roles of antioxidant enzymes in which the roles of these proteins were to remove overaccumulated ROS under stress, the main function of PPR is to inhibit the synthesis of ROS [42]. In Arabidopsis, a PPR mutant plant, in which the expression of this gene was reduced, is more sensitive to oxidative stress [42]. This is most likely caused by the lost ability on ROS generation inhibition. The increased expression of an extraordinary number of PPR genes in soybean under heat may greatly inhibit the synthesis of ROS. This could be one of the reasons for the soybean plants to cope with heat.

TABLE 4

Table 4. Differentially expressed pentatricopeptide repeat-containing protein (PPR) genes.

Accumulation of low molecular weight compounds such as proline, sugars, and polyols are often observed when plants are under stress conditions [43]. It is known that these compounds help stabilize the structure of proteins and cell membranes, buffering cellular redox, and reduce osmotic potential [44–47]. In this study, the expressions of a number of genes involved in the synthesis of these protective compounds have been found to be changed (Table 5). For example, the RING finger protein, besides the role as a ROS scavenger, was found to be also involved in the synthesis of proline, soluble carbohydrate, and osmotic regulation [48–50]. In the transgenic Arabidopsis plants with higher expression of the RING finger gene, the proline and sugar contents were increased and this is correlated to the improved resistance to stress tolerance [51]. Ornithine is a non-protein amino acid and a precursor for proline synthesis. N-acetyl-L-glutamate kinase is one of the known genes on the acetyl-glutamate pathway that is involved in the synthesis of ornithine [52]. The overexpression of the N-acetyl-L-glutamate kinase in soybean indicated that the synthesis of proline via the acetyl-glutamate pathway may be one of the mechanisms used by soybean to deal with heat [52]. Most previous genetic engineering work was focused on the manipulation of P5CS, the key enzyme for proline synthesis. However, in this study, the expression of this gene was not changed. It remains interesting to study whether the proline content can be increased in soybean plants via the overexpression of N-acetyl-L-glutamate kinase.

TABLE 5

Table 5. Differentially expressed genes related to compatible compound synthesis.

Adverse environment such as drought or heat disturbs the structure and stability of the proteins and mRNAs in plants cells. Plants have evolved mechanisms to either refold these cellular components if the “problems” were fixable or remove them if the damages were too severe. Heat shock proteins are a group of chaperone proteins that are involved in protein re-folding [53], in total, eight of heat shock proteins have been found to be differentially expressed in the heat stressed soybean (Table 6). It is fascinating to see the expression level of a 22.7 kDa class IV heat shock protein was increased by almost ten folds, indicating refolding or repairing the damaged proteins is an active process in soybean in response to heat. In addition to its negative effect on protein structure and stability, it is known that heat stress also influences some post-transcriptional regulation on mRNA which includes degradation, transport, mRNA stability, pre-mRNA processing, and localization [54]. In this study, we found three genes with functions in mRNA refolding or denaturing changed their expression levels in the heat stressed soybean plants (Table 5). For example, the expression of a DEAD-box ATP-dependent RNA helicase 53-like gene was found to be increased. The protein encoded by this gene was shown to be a mRNA chaperone that uses ATP to disrupt misfolded RNA structures and promote correct folding [55]. The altered expression of a number of genes in this category indicated that mRNA repairing or degradation is one of the mechanisms that soybean plants may use to cope with the heat stress.

TABLE 6

Table 6. Differentially expressed heat shock protein genes.

TABLE 7

Table 7. Differentially expressed genes related to E3-Uniquitin family.

Besides the “rescuing” functions derived from the heat shock proteins, a different group of proteins that participate in the degradation of damaged proteins were discovered to be differentially expressed in soybean grown under heat condition. In this study, four members of the E3 ubiquitin protein like genes were found to change their expressions, either up or down, in response to heat (Table 7). In Arabidopsis, the expression of AtCHIP, a type of E3 ubiquitin proteins, was increased in response to high temperature. However, overexpression of this gene in Arabidopsis rendered the transgenic plants to be more sensitive to temperature stress. One possibility was because high levels of AtCHIP proteins led to the turnover of some misfolded proteins before they have a chance to be folded back to functional structures [55]. The up- and down-regulation of different member of the same type of genes indicated a fine tune system exists in soybean in terms of the rate of protein turnover. It remains an interesting question to study whether the degree of damage for a protein is correlated to the expression change of their corresponding ubiquitin proteins.

Ethylene is a well-known plant hormone which not only plays crucial regulatory roles in regular plant growth and development but also serves as stress-related hormone required in various biotic and abiotic adverse growth conditions [29]. Various studies have shown that ethylene response factors (ERFs) are essential downstream components of the ethylene transduction pathway [56–59]. In Arabidopsis, AtERF1 plays a positive role in drought and heat stress tolerance by regulating stress specific gene, and by integrating jasmonic acid, ethylene and abscisic acid signals [60–62]. In transgenic Arabidopsis with overexpressed AtERF1, more compatible compounds such as proline and higher expression of heat shock protein genes were detected [60]. In tomato and tobacco, it was demonstrated that by overexpressing ERF genes, plants showed higher ethylene production and increasing tolerance to abiotic stresses. One of the functions related to ethylene under abiotic stress is its association with oxidative stress responses [63]. Ethylene might act as a signal to activate oxidative defenses during heat stress. In this study, the expression levels of 12 ethylene response factor (ERF) genes were increased (Table 8). This indicated that signaling pathway related to ethylene is involved in the response of soybean plants grown under heat. In the future, functional analysis on the identified ERFs, especially those dramatically altered one, e.g., ERF021 (78.7 fold change), in transgenic plants and the associated physiological measurements such as proline accumulation, heat shock protein gene expression, and oxidative stress status will make connections on the expression of these genes and the associated physiological outcome.

TABLE 8

Table 8. Differentially expressed genes related to ethylene biosynthesis.

CONCLUSIONS

RNA-Seq is an efficient method to study the molecular responses in soybean plants grown at heat stress condition. A total of 633 genes were found to be differentially expressed. Compared to drought stress, most of the identified genes cannot be categorized with the KEGG software. The differential involvement of antioxidant enzymes between drought and heat indicated the two stresses posed different physiological effect to the plants. Future work using genetic engineering will shed more light on what roles these identified genes may play in soybean to heat. Since this project only studied soybean plants experiencing short term heat treatment, it remains an interesting question to explore whether longer period of heat treatment will lead to stronger responses and higher number of differentially expressed genes. In addition, proteomics and metabolomics should also be used to get a systemic understanding on the physiological alterations in soybean plants experiencing heat stress because it is known that plants adjust their responses at transcriptional, translational and post-translational levels. A more complete picture on the physiological and molecular pathways plants respond to heat will provide more accurate tools for crop improvement.

SUPPLEMENTARY MATERIALS

Supplementary Figure S1: KEGG analysis of DEGs in response to heat stress (190KB)

Supplementary Table S1: KEGG-enriched pathways under heat condition (84.3KB)

DATA AVAILABILITY

The results were mean ± SE of three independent replicates. Data Availability: The transcriptomic datasets are available in NCBI with accession number SRP139514.

AUTHOR CONTRIBUTIONS

CZ, YJ, and CS conceived and designed the experiments; CX, CX, JH, MZ, SW and MI performed experiments; ZX and ZH devised and implemented the bioinformatics methodology and analysis; CX and CZ wrote the manuscript; CZ, YJ, and CS revised and finalized the manuscript.

CONFLICTS OF INTEREST

The authors declare that they have no competing interests.

FUNDING

This work was supported by funds from Purdue University as part of AgSEED Crossroads funding to support Indiana’s Agriculture and Rural Development (to C. Z.) and the USDA National Institute of Food and Agriculture (Hatch project 1007567).

ACKNOWLEDGMENTS

The authors thank the Genomics Core Facility at the Weill Cornell Medicine for their help on RNA Seq analysis and Swarup Mishra for his help on the photosynthesis measurement.