Genome and protein sequences were downloaded from the phytozome database (http://phytozome.jgi.doe.gov) and the HMMER software (http://hmmer.janelia.org) was used to build the switchgrass protein data. In addition, the Hidden Markov Model (HMM) file of PvHsp20 (PF00011) domains were downloaded from Pfam (Pfam; http://pfam.sanger.ac.uk/) (Finn et al., 2014) for the identification of Hsp20 proteins from local database (E < 0.001). All hits were confirmed by Pfam (PF00011) and NCBI Conserved Domain Search (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). The confirmed PvHsp20s were aligned using Clustal X (v 2.0) (Larkin et al., 2007) to remove the redundant sequences. The longest translated protein was picked out among PvHsp20s with alternative splicing sites, and the duplicated result was removed in phylogenetic tree analysis. The Hsp20 protein sequences in rice were referenced according to Sarkar's study (2009), while Hsp20 protein sequences of Arabidopsis were obtained from Scharf's (2001) study.

The coding sequence (CDS), exons and introns number, amino acid (aa), and chromosomal location information of switchgrass PvHsp20 proteins were retrieved from the phytozome database (http://phytozome.jgi.doe.gov). The Hsp20 proteins' molecular isoelectric point (pI) and weight (Da) were calculated by using the ExPASy program (http://web.expasy.org/compute_pi). Exon-intron display was conducted through the gene structure display server (http://gsds.cbi.pku.edu.cn). The conserved motifs among subgroups of PvHsp20 proteins were identified using the program MEME (Multiple Expectation Maximization for Motif Elicitation; version 4.11.1; http://meme-suite.org/tools/meme) with default parameters, and the maximum number of motifs to find was set to 10 for the prediction of subdomains (Bailey, 2006). The ACD sequences of ACD-containing proteins (Hsp20 and Acd proteins) were aligned via Promals3D structural alignment program (http://prodata.swmed.edu/promals3d/promals3d.php) (Pei et al., 2008). The neighbor-joining (N-J) phylogenetic tree of PvHsp20 proteins of switchgrass, Arabidopsis, and rice was built with alignments using ClustalX (bootstrap 1,000 replicates) via MEGA 5.0 (version 5.0, http://www.megasoftware.net) (Tamura et al., 2011).

The chromosomes in switchgrass were ordered to match syntenic fortail millet (Setaria italica L.) chromosome order (http://phytozome.jgi.doe.gov). We used MapInspect software (http://mapinspect.apponic.com/) (Ma et al., 2013) to generate chromosome location images to localize switchgrass PvHsp20 genes. The ratio between nonsynonymous and synonymous nucleotide substitutions (Ka/Ks) were obtained through DNAsp5 software (http://www.ub.edu/dnasp/) (Librado and Rozas, 2009).

For each of the ACD-containing proteins (Hsp20 and Acd proteins) identified in switchgrass, Unitranscript IDs were found in the PviUTs database (http://switchgrassgenomics.noble.org/) (Zhang et al., 2013). The integrated transcript sequence database was recognized through searching Unitranscript IDs in PviGEAs (http://switchgrassgenomics.noble.org/) (Zhang et al., 2013). The results from the database were graphically presented in a heatmap format using a log₂fold change after value normalization through the R Project software (http://miyoviqo.tha.im/).

Data from the ArrayExpress repository under the accession number E-MTAB-1897 (Li et al., 2013) were retrieved, for the heat-responsive transcription analysis of the PvHsp20 genes. A total of 92 ACD-containing proteins (Hsp20 and Acd proteins) retrieved from the array data were presented in a heatmap with log₂ fold change after value normalization by the R Project software (http://miyoviqo.tha.im/) (Ripley, 2001).

An interaction network of PvHsp20 proteins was built to analyze genome-wide protein-protein regulation network on the basis of orthologous rice proteins to predict the relationship between PvHsp20s and other proteins by using the rice interactions viewer (http://netbio.sjtu.edu.cn/riceppinet/search.html) (Liu et al., 2017). And then the homologs of these interactive PvHsp20 proteins in switchgrass were identified by using BLAST analysis. The interaction network of PvHsp20 proteins was drawn by Cytoscape_v3.4.0 (Smoot et al., 2011).

Switchgrass cv. Alamo seeds were sown in pots (0.2 meter diameter × 0.3 m tall) containing 1,000 g soil (pH 5.56, 1.35% organic qualitative content, 100.33 mg/kg N, 4.93 mg/kg P, and 332.25 mg/kg K). Growing in a growth chamber (Wenjiang, Sichuan, China) at 28°/20°C (day/night), the plants had photoperiod of 16 h/8 h (day/night). After germination, seedlings of switchgrass were thinned to four plants per pot. Fifty days after sowing, the potted switchgrass seedlings were exposed to a variety of stresses including ABA, drought, cold, and salt conditions as follows. For ABA treatment, the seedlings were sprayed with 100 mmol ABA for 16 days and leaves were sampled after 0, 8, and 16 days of treatment. For drought treatment, the potted seedlings were kept without watering for 28 days and at the end of drought treatment, the soil water content of drought-stressed plants was measured to be 10%. And the leaf samples were collected after 0, 14, and 28 days under drought treatment. For cold treatment, the seedlings were subjected to cold stress for 6°C for 28 days and leaves were harvested after 0, 14, and 28 days of treatment. For salinity treatment, the seedlings were watered with 250 mmol/l NaCl for 28 days, and leaves were sampled after 0, 14, and 28 days of treatment. For heat treatment, the plants were exposed to high temperature for 38°C for 28 days and leaves were collected at 0, 14, and 28 days. All materials that harvested from each treatment were frozen in liquid nitrogen immediately and stored at −80°C before for RNA isolation. With three biological replicates, all experiments were conducted three times for qRT-PCR analysis.

Total RNA kit II (Qiagen, USA) was used to isolate total RNA. RNA concentration measurement, DNaseI treatment, and cDNA synthesis were implemented as described in our previous study (Huang et al., 2014). Nine primer pairs of predicted stress-related genes were designed via Primer 5 software (Lalitha, 2004; Table S1). To confirm the primer specificity, we blasted each primer sequence to the switchgrass genome (https://phytozome.jgi.doe.gov/pz/portal.html). We also confirmed whether the nine primer pairs displayed corresponding melting curves with a single sharp peak and an electrophoresis pattern of a single amplicon with precise length. Besides, the UCE2 gene, as a reference gene, helped to calculate the expressions of nine genes (Huang et al., 2014). Moreover, the cut-off value was set to 2-fold for stress-specific expression as “down-regulated” or “up-regulated” ones.

A total of 92 proteins containing ACD were identified in switchgrass from the newly released switchgrass genome database (Panicum virgatum JGI v1.1), and detailed information about these proteins were presented in Supplementary File 1. A neighbor-joining (N-J) phylogenetic tree constructed with switchgrass, Arabidopsis, and rice proteins containing ACD domains clearly showed that 63 out of 92 switchgrass proteins were clustered with Hsp20s (or sHsps), while the rest 29 proteins with Acd proteins (Figure 1). And the 63 Hsp20s proteins were designated according to their molecular weights (Supplementary File 1). The lengths of all proteins ranged from 100 to 360aa with three of them longer than 300aa (Pavir.Ib04161.1, Pavir.J00993.1, and Pavir.J00136.1). The predicted molecular weights of the 92 ACD-containing proteins were between 11.22 kDa (Pavir.J14863.1 and Pavir.J22097.1) and 40.48 kDa (PvHsp20-40.5) with an average of 20.67 kDa. The isoelectric points (pI) of 92 ACD-containing proteins ranged from 5.04 (PvHsp20-14.6) to 11.48 (PvHsp20-15.5c), and 54 of them had pIs lower than 7.0, and 36 had pIs higher than 7.5 (noting that plant cytosol's pH is ~7.5) (Supplementary File 1). In addition, the 63 switchgrass Hsp20s were further classified into 14 subgroups in accordance with their predicated subcellular localizations: 39 Hsp20s in subgroups C (for cytosol/nucleus) I-XI, 3 in Px (peroxisome), 4 in ER (endoplasmic reticulum), 5 in P (plastid), and 12 in M (mitochondria) I-II (Figure 2).

In a simplified N-J phylogenetic tree containing only PvHsp20s, similar gene structures and motif arrangements were found among those classified in the same clade (subgroup), which consistency in turn supported the reliability of the phylogenetic classification (Figure 2). It was notable that PvHsp20s presumably located in cytosol/nucleus (CI-CX), endoplasmic reticulum (ER), and peroxisome (Px) shared strikingly different motif arrangements from those in plastid (P) and mitochondria (MI-II). And more than half of PvHsp20s (37/63) were intronless in subgroups CVIII, CIX, CXI, ER, and Px, while most proteins in subgroups CIII, CV, CVI, P, and MI-II had one or two introns.

A total of 10 types of consensus motifs were revealed among the switchgrass ACD-containing proteins, where nine types of motifs were found among 63 PvHsp20s (motif logos shown in Figure S1). Motifs-1, -2, -3, -4, and -6 were the most conserved ones (Figure 2). Matching these consensus motifs with featured sequences of ACD, we found that the Consensus Region I (CRI) of ACD was composed of motifs 2, 4, and 6, while CRII was composed of motifs 1 and 3. The presence of all of these five motifs composing the ACD were only found in one half (31 out of 63) of PvHsp20s in subgroups of CI, CII, CVIII, CIX, CX, and ER, while the rest PvHsp20s lacked at least one of these five motifs, suggesting the diversification with the ACD across PvHsp20s, which was proposed as one key domain for protein-protein interaction.

Multiple alignments of the conserved ACD domains among PvHsp20s were conducted as shown in Figure 3 to illustrate the relationship between consensus motifs and predicated secondary structures of ACDs. The CRI of ACD were composed of motifs 2, 4, and 6 featured with four β-sheets (β2, β3, β4, β5), and CRII were composed of motifs 1 and 3 featured with three β-sheets (β6, β7, β8, and β9). Among the 63 PvHsp20s, only PvHsp20-11.6 lacked β7, PvHsp20-15.0b lost β5, and HSP17.5a missed β8 and β9 (Figure 3). Additionally, a total of 10 PvHsp20s (15.9%) lacked β6. Because β6-strand is significant for dimer formation and oligomerization by strand swapping (Van et al., 2001), the process of dimer formation or oligomerization for these PvHsp20s without β6-strand might be influenced. These PvHsp20s lacking of β6 might depend on other courses to remedy this omission for protein-protein interaction (Waters, 2013).

Heat shock protein (HSP) protein families play a significant role in heat stress tolerance. With reference to functional annotated Hsp20s (Table S2), we identified 28 corresponding orthologous PvHsp20s that were potentially involved in the regulation of abiotic stress tolerance in nine functional groups (Figures 4A–I).

Furthermore, the expression profiles of PvHsp20 and Acd genes from a switchgrass Affymetrix array (Li et al., 2013) were analyzed to discover those transcriptionally responsive to heat stress. A total of 55 PvHsps (Figure 5) and 22 Acds (Figure S2) were discovered in the switchgrass GeneChip (Supplementary Files 2, 3). Among the 55 PvHsps, 48 were up-regulated, and 24 of them were with more than 1.5-fold change and four with more than 5-fold change after heat stress. As for the Acd genes, nearly half of them were down-regulated (9/22, 40.9%), only three were up-regulated, and the rest 11 genes did not display an obvious change (below 1.2-fold change after heat stress) after heat stress. These results were consistent with previous report that most PvHsp20s were transcriptionally up-regulated and some Acd genes were down-regulated by heat stress (Sarkar N. K., 2009).

Among the 24 up-regulated PvHsp20s with more than 1.5-fold change after heat stress, seven of them (PvHsp20-17.4b, -23.1a, -23.1b, -19.5b, -17.8b, -24.2, and -23.7) were orthologous to previously functional annotated Hsp20s in model plant species (Table S2), further suggesting that these PvHsp20s might be involved in plant heat tolerance.

In the absence of stress, differential expression of Hsp20s during development stages and/or in different tissues has also been recorded before (Vierling, 1991). Therefore, PvHsp20s' expression patterns in 21 different tissues and developmental stages were analyzed using switchgrass Gene Expression Atlas (PviGEA). According to the analysis, a total of 39 PvHsp20s had differential expression patterns in different organ/tissues. (Figure 6). Specifically, 38 genes displayed relatively higher expression levels in reproductive organs (inflorescence and seeds at different developmental stages) (Figure 6). These accumulations of PvHsp20s might contribute to the thermotolerance in reproductive organs (florets and seeds), and indicated that PvHsp20s' significance in normal plant growth and development even under optimal growth condition. As for the Acd genes, none of them was found with tissue-specific expression pattern (Figure S3).

To our knowledge, there was no report concerning Hsp20s' involvement in secondary cell wall strengthening and/or lignification. However, the gene expression altas showed that one PvHsp20 gene (PvHsp20-19.2) was highly expressed in lignified tissues (e.g., crown, roots, node, internode, and inflorescence branches) (Figure 6). Concerning the importance of lignin content in this lignocellulosic biomass grass, it would be interesting to further understand the function of PvHsp20-19.2.

Based on in silico data analyses described above, expression profiles of nine PvHsp20s that were predicted as stress-related genes and/or orthologues of functional-annotated Hsp20 genes, were further studied under drought, ABA, salt, cold, and heat treatments using qRT-PCR (Figure 7; Supplementary File 4). Setting the cut-off value at 2-fold change, eight out of nine selected PvHsp20 genes were expressed at significantly higher or lower levels when treated with at least one of these abiotic stresses. Among these PvHsp20s, only three genes (PvHsp20-16.1, -17.4, and -21.8b) were transcriptionally up-regulated by heat stress, and another three (PvHsp20-17.4b, -19.2, and -19.5b) down-regulated by heat stress, respectively. On the other hand, all of the tested PvHsp20s were responsive to drought treatment, among which five were significantly up-regulated by under severely drought stress (drought treatment for 28 days), one was up-regulated under moderate drought stress (drought treatment for 14 days), and three were significantly down-regulated under moderate drought condition. When exposed to salt treatment, five genes containing PvHsp20-16.1, -17.4a, -17.5b, -19.5b, and -21.8b were significantly up-regulated after 14 and 28 d, while only one gene (PvHsp20-21.0b) were significantly depressed. Cold treatment significantly induced expression of six out of nine tested genes except PvHsp20-16.7b, -17.4b, and -21.0b. In particular, PvHsp20-19.2 and -19.5b increased over 20 fold change after 14 and 28 d of cold treatment. Yet, the up-regulation of PvHsp20 genes to these abiotic stresses seemed to be independent on ABA that nearly all of these were down-regulated in response to ABA treatment, except PvHsp20-21.0b.

PvHsp20s function as chaperones. Therefore, an interaction network of PvHsp20 proteins was built on the basis of orthologous rice proteins to predict the relationship between PvHsp20s and other proteins in switchgrass (Figure 8; Supplementary File 5). Nearly all abiotic stress-related PvHsp20s as described above (84.6%, 11/13) were built in this network. A total of 109 high confidence interactive relationships (score > 0.85) and 68 interactive proteins were discovered, including 13 PvHsp20s and 55 other switchgrass proteins. The interaction network of PvHsp20 proteins displayed a complex functional relationship. All PvHsp20 proteins directly or indirectly interacted with each other, and according to the results from Table S2 and affymetrix array (Figure 5). PvHsp20-16.9d was predicted to directly interact with the largest number (35) of proteins in switchgrass, followed by PvHsp20-17.8d (27), PvHsp20-15.5b (22), and PvHsp20-19.1 (13), indicating their core status of stress regulation in the protein interaction. In addition, PvHsp20-15.5b, -19.1, -17.8d, and -16.9d all interacted with one drought inducible protein (Pavir.J33423.1) encoding a DnaJ superfamily chaperone (Ye et al., 2016).

Allotetraploid switchgrass has two subgenomes designated as A and B (Okada, 2010). In this study, chromosomal localizations of 41 PvHsp20s were further analyzed. As shown in Figure 9, these genes were unevenly distributed on 14 chromosomes of five homeologous pairs (Chr01a/b, Chr04a/b, Chr05a/b, Chr06a/b, and Chr09a/b) and two nonhomologous chromosomes (Chr02b, Chr03a, Chr07a, and Chr08b). For example, there was only one gene (14.63%) located on chromosomes 02b, 03a, 06a/06b, 07a, and 08b; while there were seven genes (17.07%) on chromosome 05a. A total of 6 pairs of paralogous ACD-containing genes (bootstrap value >95 in the phylogenetic tree) with defined chromosomal locations were linked with red line in Figure 9, and all of them were in homeologous chromosomes (Table 1). Tandem gene duplication, defined as paralogous genes physically linked in tandem with less than five genes in-between, was not found in this study, suggesting that these paralogous ACD-containing genes were all duplicated due to the allotetraploidy event of switchgrass. During the time course of evolution, favorable mutations leading to species divergences were usually fixed (called “diversifying selection”), while those causing disadvantages were eliminated (called “purifying selection”). Comparison of synonymous and nonsynonymous substitution rates between the six pairs of PvHSP20 paralogous genes suggested that 83.33% (5 out of 6) gene pairs were under purifying selection, while only one pair (1/6, 16.67%) was under diversifying selection (Table 1).

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.