{'fraction_non_alphanumeric': 0.03632246804094533, 'fraction_numerical': 0.044954950705221944, 'mean_word_length': 4.870949875380781, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Expression of Heat Shock Proteins by Heat Stress in Soybean

Heat stress is one of the factors disturb productivity and growth of plants. Many genes including heat shock protein (HSP), heat shock transcription factors (HSF) and chaperones, were identified and characterized in many plants to play role in increased tolerance to abiotic stress. To reveal responsive gene to heat stress, we performed RNA-seq using two Korean soybean varieties under heat stress and normal conditions. The transcripts were analyzed, and we obtained 2,458 genes including 46 co-up regulation and 55 co-down regulated genes in both soybean varieties. We also revealed HSPs, HSFs and chaperones in the differentially expressed genes using BLAST and Pfam analyzation and verified expression changes under heat stress. Finally, we find 68 genes involved in HSP, HSF, chaperones in heat responsive genes associated increasing heat tolerance. As a result, relatively small HSP families were up regulated and continuously expressed in long period heat stress. On the other hand, large molecule HSPs, HSFs and chaperonin did not response to long heat stress. The expression profiling and characterization provide invaluable information to understand heat tolerance of soybean.


INTRODUCTION
Global climate change is becoming a serious condition to plants. Among abiotic stress, heat stress is one of the factors disturb productivity and growth of plants. Plants could not move their position to avoid stresses, stress response system with molecular level is important to growth and productivity in plants. Among many stress response systems, heat shock proteins (HSPs), heat shock transcription factors (HSFs) and chaperone are essential for maintenance and restoration of protein structure and stabilized condition of plants under heat stress (Boston et al. 1996;Nakamoto and Vigh 2006;Xu et al. 2012). Since the HSPs expression observed under heat stress (Lindquist and Craig 1988), the heat stress response pathway has been continually studied in plant (Queitsch et al. 2000;Kotak et al. 2007;Larkindale and Vierling 2007;Yamada et al. 2007). Heat stress has make mis-folded proteins and cause to losing of biological activity in plant. These damaged proteins are interfering to normal cell functioning. Chaperones and HSPs have functions which restoration of mis-folding proteins and minimizing damages for maintenance plant homeostasis (Boston et al. 1996;Bukau et al. 2006;Pratt et al. 2010). Also, these stress response systems are regulated by HSFs which has DNA binding sites (Schöffl et al. 1998;Nover et al. 2001;Hsu et al. 2010).
Plant HSPs are classified into some types according to molecular size (Sarkar et al. 2009), small HSPs, HSP20, HSP60, HSP70, HSP90 and HSP100. These families also have been characterized in a few plants (Krishna and Gloor 2001;Hu et al. 2009;Lopes-Caitar et al. 2013). Some studies had tried to reveal a function of HSPs, HSFs and chaperones from Zea mays, rice, Arabiodopsis and soybean (Queitsch et al. 2000;Rhoads et al. 2005;Liu et al. 2009;Qi et al. 2010;Sun et al. 2012) and expression gene profiles of HSPs and HSFs were identified using microarray analysis and quantitative real time PCR (Hu et al. 2009;Lopes-Caitar et al. 2013).
Since the development of next generation sequencing method, gene and transcript analysis has made possible with quickly and cost effectively. Because of the low cost and high capacity of data, RNA-sequencing (RNA-seq) has been using to measuring differentially expressed genes (DEGs) in response to stress (Marioni et al. 2008;Wang et al. 2009;Anders and Huber 2010;Rapaport et al. 2013). It was also used to reveal differential expressed transcriptome in rice under cold conditions (Shen et al. 2014), HD-Zip transcription factors in soybean under salt stress (Belamkar et al. 2014) and regulated genes in Brassica juncea under drought stress and high temperature (Bhardwaj et al. 2015).
Soybean is one of important crop plant in worldwide. In soybean, heat stress significantly decreased plant growth and yield. Recently, mechanism of HSP20 and HSP90 families were studied to understand of heat response system in soybean (Lopes-Caitar et al. 2013;Xu et al. 2013). But there were not enough attempts to reveal the profiling of DEGs involved in HSP families under heat stress in soybean.
In this study, gene expression pattern was analyzed from RNA of soybean leaves under normal and heat stress conditions. The plants were treated for 12 hours under heat stress to focus on verifying referential genes to increase heat tolerance. We used RNA-seq as a tool to obtain differential expression pattern of genes under heat stress. The DEGs were analyzed by Blast to identify their functions, and compared to gene lists of HSPs, HSFs and chaperons were obtained from Blast searches on Phytozome. This study not only provides heat responsive genes but also improves knowledge for HSPs, HSFs and chaperons.

Plant materials and stress conditions
Seeds of soybean accessions Daewon (IT 212859) and Hwangkeum (IT 142807) were obtained from National Institute of Crop Science, RDA, Suwon, Korea. These seeds were planted directly into pots containing Sunshine mix 5 (Sungro Horiticulture, Quicy, MI) and were grown until V5 stage in greenhouse. After V5 stage was reached, the plants were moved to heat stress conditions. The plants were exposed to heat stress (38 ± 2°C) for 12 hours and then leaves were immediately collected in liquid nitrogen gas and transferred to deep freezers (-80°C).

RNA isolation and RNA-seq
For all samples, total RNA was prepared from leaf tissues under each condition using an RNeasy Plant Mini Kit (Qiagen, Germantown, MD, USA). The RNA quality was checked for integrity before performing the RNA sequencing process using a 2100 Bioanalyzer RNA Nanochip (Agilent Technologies, USA). After a quality check, RNA sequencing using Illumina HiSeq platform was performed by Theragen, Inc. (Korea). Raw sequence reads have been submitted to the NCBI sequence read archive (SRA) database under accession number SRA274425.

RNA-seq data processing and annotation
Raw reads from Illumina HiSeq 2000 sequencing platforms were imported to CLC Genomics Workbench 7.5.1 (CLC Bio) and read statistics were assessed using sequencing data quality control, followed by read trimming for quality. The raw reads were trimmed with a quality score limit of 0.01. The high-quality reads were mapped to the soybean genome sequence of Williams 82 from Phytozome version 10 (Gmax_275_wm82.a2.v1) (http:// www.phytozome.net) used as the reference sequences in CLC Genomics Workbench. The RNA-seq analysis was carried out gene expression and transcript expression using default parameters.
The gene expression data of either in normal or heat stressed plants was determined using reads per kilobase of transcript per million mapped reads (RPKM) used as a measure (Mortazavi et al. 2008). The expression values were normalized using scaling method to comparable (Robinson and Oshlack 2010). After normalization of the RPKM values, the expression fold change was calculated based on normalized RPKM log2 transformed expression values between each condition. We used Kal's method to test the significance of the expression fold change (Kal et al. 1999). To select differentially expressed genes, we set threshold as P-value ＜ 0.05, mapped reads ＞ 5 and fold  (Conesa et al. 2005) was used to annotate each consensus sequence of the assemblies. Each consensus sequence was analyzed using BLAST to nonredundant (nr) database of NCBI using BLASTX with e-value of 1.0E-3 and default parameters. The consensus sequences contained Hidden Markov Model (HMM) domains were searched in Pfam databases (version 27). The results of Blast2GO were applied using the hyper-geometric tool BiNGO as cytoscapes app (Maere et al. 2005).

HSPs, HSFs, chaperones analysis
We got a total number of 392 genes from BLAST results of Phytozome (version 10). The list of BLAST results contained 276 HSP genes, 52 HSF genes and 64 chaperone genes. The 392 genes were compared with expression values from RNA-seq data and selected expressed genes under heat stress. The consensus sequences of selected genes were conducted to phylogenetic analyses using neighbor joining distance algorithm with 1000 bootstrap.

Validation of differentially expressed genes from RNA-Seq using quantitative real time RT-PCR analysis (qRT-PCR)
The extracted RNA was reverse transcribed to cDNA using PrimeScript 1 st strand cDNA synthesis kit (Takara Bio, Otsu, Japan) and random primers according to the manufacturer's procedure. The cDNA was performed on Rotor-gene Q real-time PCR system with SYBR Green PCR Kit (Qiagen, Germantown, MD, USA) and transcript specific primers. The cDNAs have been checked expression of heat stress marker (GMHsp90C2.1) which can validate for the plants applied heat stress (Xu et al. 2013

Analysis of RNA sequencing data
To compare expression profiles of genes, we obtained two pairs of RNAs from control and heat-treated plant in Daewon and Hwangkeum. Total number of raw reads ranged from 25,678,226 to 27,543,208 in Daewon and Hwangkeum has shown reads range from 51,460,028 to 51,764,174 (Table 1). After trimming and filtering, approximately 98.1 million reads aligned to the reference genome (Gmax_275_wm82.a2.v1). Of the total paired reads, average over of 60% reads were mapped to reference genome. Among mapped reads, 92.2-95.7% reads were located to exons. Some reads were located to introns (1.0-1.3%) and intergenic regions (3.3-6.5%) (Supplementary Table S1). The genes with mapping reads were used for differentially expressed genes (DEGs) analysis.

Differentially expressed genes (DEGs) analysis
Totally, the abundance of 40,474 genes and 41,845 genes was shown as a result of RNA-seq under heat stress in Daewon and Hwangkeum, respectively. To obtain heat-regulated genes, we selected differentially expressed genes by comparing normalized RPKM log2 transformed expression values. After the comparing, 1,292 and 1,166  genes were obtained between heat-stressed plant and control plants with thresholds in Daewon and Hwangkeum, respectively (Supplementary Tables S2 and S3). A total number of 2,458 genes were shown differential expression pattern. A number of 347 genes were up regulated by heat stress in Daewon and 945 genes in down regulated. The down regulated genes were found 644 genes, while 522 genes were up regulated under heat stress in Hwangkeum. Among them, the 46 genes had been found as co-up regulated genes under heat stress conditions and the 55 genes were shown co-down regulated genes (Fig. 1). A total number of 51 genes were identified inversed regulation between two varieties.

Gene ontology (GO) enrichment analyzation
We selected 1,292 genes and 1,166 genes in Daewon and Hwangkeum by DEGs analysis, which include co-up and co-down regulated genes. The DEGs were assigned to functional categories using Blast2GO. The DEGs were  Table  S4). As remarkable results, many genes were classified to transcription factor activity (GO:0003700, 91 sequences), transcription regulator activity (GO:0030528, 92 sequences) and sequence-specific DNA binding (GO:0030246, 58 sequences). These GO terms were reported helpful to improve abiotic stress tolerance including HSP and HSF

HSPs, HSFs and chaperones
We obtained gene list of HSPs, HSFs and chaperones from Blast searches in Phytozome. The list contained 276 HSP genes, 52 HSFs and 64 chaperone genes from soybean. Among them, only 68 genes were shown within DEGs ( Table 2). The 5 genes had shown as co-up regulated genes, contained 4 HSP genes (HSP20/alpha crystallin family) and 1 AAA domain (Cdc48 sub-family) have reported play a role as like HSP (Supplementary Table S5). By contrast, the only 1 gene was shown in list of co-down regulated gene. The expression patterns of HSP, HSF and chaperonin are different for cultivars. Daewon has shown down-regulation of HSP families and up-regulation of chaperonin under heat condition. But relatively more HSPs are shown up regulation to response heat stress in Hwangkeum. Generally, HSP70 families were down regulated and HSP20 families were up regulated in both cultivar under heat stress (Supplementary Tables S5 and  S6). A large part of the co-up regulated genes is HSP20, in comparison, the co-down and down list each variety were consisted of HSP70 and HSP90 ( Fig. 3 and Supplementary  Table S5). Likewise, HSP20 families were shown in list of up-regulated genes from each variety, and large molecular HSP families (HSP70 and HSP90) were shown downregulation. Almost the HSFs were not found in up regulated genes, but HSFs were found in down regulation list of each variety.

Validation of differentially expressed genes by quantitative real time RT-PCR
To verify reliability of the expression profiles from RNAseq data, qRT-PCR analysis was performed for DEGs (Fig.  4). We have checked the expression of heat stress marker (GMHsp90C2.1) with three biological replicates. Expression of GMHsp90C2.1 which has been used as heat stress marker was shown up-regulation in both cultivars. As a result, it was confirmed that the plant was treated with heat stress. RNA-seq results showed that almost of the large molecule HSPs were down-regulated, and the same results were obtained with qRT-PCR. Co-down regulation gene (Glyma.18G289100) and four genes (Glyma.08G025700, Glyma.18G074100, Glyma.13G224000 and Glyma. 15G088000), which showed down regulation in each cultivar showed similar pattern in qRT-PCR. The other hand, HSP20 gene family increased in both cultivars under heat stress. In both cultivars, 5 genes (Glyma.04G054400, Glyma. 06G054700, Glyma.07G200200, Glyma.08G242100 and Glyma.18G094600) were expressed as co-up regulation, which could be confirmed by qRT-PCR. The other genes (Glyma.04G229800, Glyma.09G131500 and Glyma. 17G072400) that were up-regulated in each variety were not distinguished by DEGs because of the high cut-off of this experiment, but they showed a sufficiently high expression level. HSF genes also were confirmed by qRT-PCR. The HSF genes were shown similar results to expectation from RNA-seq while showed inconsistent expression pattern with distance of phylogenetic tree. Two chaperonin genes (Glyma.15G250500 and Glyma.17G047200) were shown similar expression pattern. Almost chaperonin genes were up regulated in Daewon while the change of expression pattern could not be observed in Hwangkeum. Generally, results of qRT-PCR nearly matched expectation from RNA-seq data. Expression pattern of DEGs which other than HSP genes were also like expected values.

DISCUSSION
The heat stress has negative effect on the physiology and metabolism of plant cells. In this study, we obtained DEGs on heat treated Korean soybeans involved in heat tolerance to verify response genes in soybean. Also, HSPs, HSFs and chaperone genes which associated with heat stress tolerance were found in DEGs. The DEGs were performed Blast and Pfam analysis for understand heat response system in soybeans.
The primary subject of this study was to observe the differentially expressed genes under heat stress. Totally, 2,458 genes were selected as DEGs in both cultivars (Fig.  1, Supplementary Tables S2 and S3) and the genes were classified according to gene ontology (GO) (Fig. 2,  Supplementary Table S4). The smaller numbers of genes which has same expression patterns as like co-up (46 genes) and co-down (55 genes) were detected than expectations. These results were due to the following reasons: First, the RNA-seq data were analyzed with highly cut-off to obtain more accurate information. For example, the HSP20 families were up regulated in both cultivars. In the alpha crystallin family, which is part of hsp20, it consists of a total of six genes (Glyma.06G054700, Glyma. 20G015900, Glyma.04G054400, Glyma.14G100000, Glyma.14G099900 and Glyma.17G224900). Due to the high cut-off of the experiment, some genes did not find to be DEGs, but expression patterns are all up-regulation. Second, genetic properties of soybean cultivars. Expansion of the HSP gene family has generally occurred through and whole-genome duplications events. We also found this phenomenon in the phylogenetic tree of DEGs. The combined phylogenetic tree revealed that the orthologous relationships between HSP families are contained in paralogous pairs of each cultivars. In conclusion, although one gene does not show the same expression pattern in both cultivars, gene families with similar roles have similar expression patterns. For example, Glyma.02G244300 and Glyma.14g214800 were up-regulated in only one cultivar. Both genes are expected to play a similar role because of same gene family, but only one gene is expressed due to their genetic characteristics. This phenomenon can be regarded a difference due to genetic diversity as a characteristic of a plant having a variety of substitution genes on one function.
To verify expression pattern of HSPs, HSFs and chaperones, we obtained list from Blast in Phytozome database. We compared the list and DEGs from RNA-seq data, and found out 68 heat responsive genes included 56 HSPs, 7 HSFs and 5 chaperones. Of the total, HSP20 families had a large part in co-up regulated genes. The HSP20 families had known play a role as protect protein denaturation in eukaryotic and prokaryotic cells (Feder and Hofmann 1999;Cashikar et al. 2005). We found HSP20 families in co-up regulated genes (Glyma.06G054700, Glyma.04G54400, Glyma.07G200200, Glyma.18G094600 and Glyma.08G242100) which were similar with the previous study (Lopes-Caitar et al. 2013). The other HSP20 families were up regulated in both cultivars, but were not passed cut-off in Daewon. Some of HSP90 genes include GmHsp90A1 (Glyma.09g131500) and GMHsp90C1.2 (Glyma.02g124500) were reported as rapidly accumulate and be maintained for a long period (Xu et al. 2013). These genes were found in up-regulation genes of Hwangkeum, but did not found to be DEGs in Daewon. As above described, this the gene up regulated in Daewon, but did not passed our high cut-off. The other HSP90 and HSP70 families were reported that rapidly accumulated and decreased expression level with time under stress (Swindell et al. 2007). Likewise, the other HSP90 families were shown down regulation in both varieties. Other large molecular HSP families as like HSP70 (Glyma.18G289100, Glyma.11G140500, Glyma.08G025700, Glyma.05G219600, Glyma.08G025900, Glyma.12G064000, Glyma.15G088000 and Glyma.13G224000) and HSP90 (Glyma.18G074100, Glyma.08G332900 and Glyma.17G258700) were down regulated by heat stress. It is seemed to that these genes had responded at short period under heat stress and decreased following time.
The HSF genes also were response to short term heat stress . Some HSF genes were showed highly expression under heat stress until 6 hours and the expression decreased at 12 hours. In our results, very few HSF genes were found. It is suggested that HSFs does not significantly affect long-term stress over 12 hours.
Gene ontology analysis was used to verify domain of co-up regulated genes as like candidate of gene related heat responses. Some genes (Glyma.03G130400, Glyma. 03G130600, Glyma.08G215500, Glyma.12G178500, Glyma.01G198100, Glyma.13G251300 and Glyma. 19G222000) with DNA binding domains (HLH, and Zinc binding) were also found in Pfam analysis. These domains were predicted to play a role as transcription factor as like HSF genes (Downs and Heckathorn 1998;Baniwal et al. 2004). Chytochrome P450 domains found in up regulations in both cultivars. This domain was known that induced by abiotic stress and contained the recognition site of TGA-box, W-box, MTB and MYC for transcription factors (Narusaka et al. 2004). Some of the other domains (FAD binding and NADH (P) binding) were seemed to related in reactive oxygen species (Supplementary Table S6) (Marino et al. 2012;Kikusato and Toyomizu 2013). These domains have a role control expression of genes as like HSFs, and could lead to improve heat stress tolerance of plants.
In conclusion, we performed RNA-seq of two Korean soybeans to understand the gene expression and responsive system under heat stress. The expression patterns and characters were determined using RNA-seq and bioinformatics tools for verify fundamental genes to response to heat stress. Based on the DEGs and Pfam analysis, the genes were inferred which the primary genes to improve heat tolerance under heat stress. Especially, HSP20 and HSP90 families could played an important role and improved tolerance in long period heat stress and maintain their homeostasis. Among the DEGs expected to play an important role, some genes did not express in same pattern in both cultivars. However, we have discovered a phenomenon that expression of similar genes to be maintain the same function. These results suggest that the heat stress tolerance of a plant is related to the expression pattern of the HSP gene families rather than to a specific gene. These results can be further analyzed to reveal gene related heat tolerance and understand responded mechanism to heat stress.