Genome-Wide Analysis and Expression Profiling of the Heat Shock Factor Gene Family in Phyllostachys edulis during Development and in Response to Abiotic Stresses

Heat shock transcription factors (Hsfs) play crucial roles in regulating plant responses to heat and other stresses, as well as in plant development. As the largest monopodial bamboo species in the world, how to adapt to various stresses under the background of global climate change is very important for the sustainable development of bamboo forest. However, our understanding of the function of Hsfs in moso bamboo (Phyllostachys edulis) is limited. In this study, a total of 22 non-redundant Hsf genes were identified in the moso bamboo genome. Structural characteristics and phylogenetic analysis revealed that members of the PheHsf family can be clustered into three classes (A, B and C). Furthermore, PheHsfs promoters contained a number of stress-, hormoneand development-related cis-acting elements. Transcriptome analysis indicated that most PheHsfs participate in rapid shoot growth and flower development in moso bamboo. Moreover, the expression patterns of all 12 members of class A were analyzed under various stresses (heat, drought, salt and cold treatment) through Figurereal-time quantitative polymerase chain reaction (qRT-PCR). Within the class A PheHsf members, PheHsfA1a was expressed mainly during moso bamboo development. Expression of four PheHsfA4s and one PheHsfA2 (PheHsfA4a-1, PheHsfA4a-2, PheHsfA4d-1, PheHsfA4d-2, and PheHsfA2a-2) was up-regulated in response to various stresses. PheHsfA2a-2, PheHsfA4d-1 and PheHsfA4d-2 were strongly induced respectively by heat, drought and NaCl stress. Through co-expression analysis we found that two hub genes PheHsfA4a-2 and PheHsfA4a-1 were involved in a complex protein interaction network. Based on the prediction of protein interaction networks, five PheHsfAs (PheHsfA4a-1, PheHsfA4a-2, PheHsfA4d-1, PheHsfA4d-2, and PheHsfA2a-2) were predicted to play an important role in flower and shoot development and abiotic stress response of moso bamboo. This study provides an overview of the complexity of the PheHsf gene family and a basis for analyzing the functions of PheHsf genes of interest.


Introduction
Moso bamboo (Phyllostachys edulis (Carrière) J. Houzeau, synonym Phyllostachys heterocycla (Carrière) is a large woody bamboo of high ecological, economical and cultural value in Asia.Under suitable spring conditions, its shoot can grow from 0 to 20 m in 45-60 days [1].Moso bamboo forest covers an area of 3.87 million hm 2 , accounting for up to 70% of the Chinese bamboo forest area [2,3].Because of its rapid growth and highly lignified culms, the annual economic value of moso bamboo production, including Forests 2019, 10, 100 2 of 17 timber and wood production, reaches 184 billion dollars [4].Moreover, carbon sequestration in moso bamboo is two to four times greater than that of Chinese fir, making it an important global non-timber forest resource [5].The growth of bamboo is dependent on natural precipitation and is vulnerable to high temperature and drought.Liu et al. [6] has shown that temperatures >40 • C and drought for >10 days during August result in severe losses in moso bamboo forests.Drought during spring can reduce moso bamboo shoot growth, yield and quality.From July to September, high temperatures and drought affect the sprouting phase of bamboo.These stresses affect the yield and quality of winter shoots, as well as new bamboo yield in the following year and the yield of wood during subsequent years [7].Climate change has also been associated with more frequent high temperatures and drought, which in turn reduce the ecological and economical value of moso bamboo.Therefore, it is essential to elucidate the molecular mechanisms involved in heat and drought stress responses to improve stress tolerance in moso bamboo.
To survive high temperature and other stresses, plants have evolved a series of defense strategies [8,9].Heat shock proteins (HSPs) act as molecular chaperones that protect cells against heat and other stress damage by preventing protein aggregation [10,11].As the terminal components of the stress signal transduction chain, heat shock stress transcription factors (Hsfs) bind to the promoter regions of HSP genes to regulate transcription in response to stress [12,13], particularly heat stress [14].Hsfs also contain a highly conserved DNA-binding domain (DBD) at their N-terminal and an oligomerization domain (OD or HR-A/B region) composed of two hydrophobic heptad repeats.Based on the amino acids of the HR-A/B region, plant Hsfs are grouped into three main classes (HsfA, HsfB and HsfC) [15].Certain Hsfs contain a nuclear localization signal (NLS), a nuclear export signal (NES) and an activator motif (AHA).The AHA motif located at the C-terminus in class A Hsfs confers transcriptional activation.Moreover, a repressor domain (RD) that contains the tetrapeptide LFGV occurs at the C-terminal of class B Hsfs.
Recent studies have revealed that plant Hsfs plays important roles in generating responses to heat and other stimuli, as well as in organ development [16].Class A HsfA1a is regarded as a master regulator and has a unique role in eliciting heat stress responses in tomato (Solanum lycopersicum L.) [17].HsfA2 is functionally similar to HsfA1 in regulating thermotolerance, as well as serving as a key regulator in osmotic and oxidative stress responses [18][19][20].The expression of HsfA3 is induced in Arabidopsis by heat and drought stress, indicating that HsfA3 might play a role in drought and heat stress signaling [21,22].The thermotolerance of Arabidopsis with hsfA3 T-DNA insertion mutants was decreased [23].Moreover, the ectopic overexpression of SIHsfA3 increases thermotolerance and salt hypersensitivity in germination in transgenic Arabidopsis [24].The Arabidopsis mutant athsfa4a was more sensitive to dehydration.Furthermore, desiccation tolerance was rescued in athsfa4a/BnHSFA4a seeds to similar levels compared with those of Col-0 [25].Transgenic chrysanthemum overexpressing CmHSFA4 displayed enhanced salinity tolerance partly due to enhanced Na + /K + ion and ROS homeostasis [26].Interestingly, in rice, wheat and Sedum Alfredii, HsfA4a possibly confers cadmium tolerance [27,28].In addition, AtHsfA9 plays a role in embryonic development and seed maturation in Arabidopsis [29].In group B, the majority of HsfBs act as repressors due to the RD region [15].However, AtHsfB1 act as repressors of the heat shock response under non-heat-stress conditions, but act as positive regulators of heat shock proteins under heat-stress conditions [30].In group C, FaHsfC1b from Festuca arundinacea confers heat tolerance in Arabidopsis [31].
Although several Hsf family genes in Arabidopsis and other plant species have been characterized, functional analysis of those in moso bamboo has been limited.The completion of the draft genome sequence of moso bamboo has greatly facilitated the identification of Hsf family at the whole-genome level [32].
In this study, we describe the genome-wide identification and analysis of the PheHsf family of moso bamboo for the first time.In addition, expression patterns of PheHsf genes during development, as well as in response to various abiotic stresses, were also investigated.Our results will provide a foundation and valuable information for future functional analysis of PheHsfs.

Database Searches for Hsf Genes in Moso Bamboo and Analyses of Physicochemical Characteristics
For more accurate identification of Hsf genes in Moso bamboo, multiple database searches were performed according Hou et al. [33].First, Hsf mRNA sequences of Oryza sativa and Brachypodium distachyon, obtained from NCBI Nucleotide database (https://www.ncbi.nlm.nih.gov/) as query sequences to blast against the moso bamboo genome database.For the filtration step of the blast process, the Hsf genes were obtained by blast bamboo transcriptome using Hsf mRNA from some other species identified as query with loose e-value of <0.00001.Next, the protein sequences of putative genes obtained from the first step were blast against the NCBI non-redundant protein database with e-value of <0.0000000001 by Blast2GO to confirm the identification.The putative genes described as Hsf proteins or proteins belonging to Hsf family were kept, and genes described as other family proteins were abandoned.Then the HSF domain (PF00447) of the HSF family were researched in these putative Hsf proteins to reconfirm our data using the Pfam database (http://pfam.sanger.ac.uk/) and Conserved Domains Database (http://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml) with e-values <0.001.Only the gene models containing HSF domain were considered to belong to the Hsf family.Finally, the selected Hsfs were further screened using the full-length non-chimeric (FLNC) reads (http://www.forestrylab.org/db/PhePacBio/)[34].The moso bamboo Hsf genome sequences, coding sequences, protein sequences, and putative novel or mis-annotated Hsf genes were obtained from moso bamboo genome database.The amino acid sequences of Arabidopsis, rice (Oryza sativa), and B. distachyon Hsf proteins were downloaded from Plant Transcription Factor Database v4.0 (PlantTFDB 4.0, http://planttfdb.cbi.pku.edu.cn/) and Heatster (http://www.cibiv.at/services/Hsf/)[35,36].Molecular weights and theoretical isoelectric point (pI) were determined using ExPASY (http://web.expasy.org/compute_pi/).CELLO v2.5 Server (http://cello.life.nctu.edu.tw/) was used to predict the protein subcellular locations for candidate PheHsfs [37]

Phylogenetic Analysis
Arabidopsis and rice Hsf gene datasets [38] were used to classify the moso bamboo Hsf genes and predict their functional roles.Multiple sequence alignments of full-length Hsf proteins from Arabidopsis, rice, B. distachyon, and moso bamboo were performed using ClustalX 1.83 (http:// www.clustal.org/)and two online programs (Clustal Omega and MUSCLE) [27,[39][40][41]].An unrooted neighbor-joining (NJ) phylogenetic tree was constructed using MEGA7.0 with 1000 bootstrap replicates.Another phylogenetic tree with only the Moso bamboo Hsfs was also constructed using the amino acid sequences according to the same method.

Cis-Regulatory Element Analysis of PheHsf Genes
The 1500-bp sequence upstream from the initiation codon of each PheHsf gene was obtained from the moso bamboo genome database.These sequences were used to identify cis-acting regulatory elements with the online program Plant CARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

Plant Material
After the surface was sterilized with 1% formaldehyde, the moso bamboo seeds were germinated in Petri dishes (12-cm diameter) lined with filter paper and containing 10 mL of sterile water.After 4 days, the germinated seedlings were planted in vermiculite and watered with 1/2 Hoagland's nutrient medium weekly (all plants were grown in 16 h day/8 h night at 22 • C).Two-month-old seedlings were used for abiotic stress treatments.According to Cheng et al. [44], drought and stress was conducted by incubated seedlings with 20% (m/v) PEG6000 and 200 mM NaCl, respectively.Heat stress and low temperature treatments were respectively created by placing seedlings in a 42 • C and 4 • C lighted growth chamber according to Liu et al. [45].The control seedlings were grown without any stress treatment.The second or third mature leaves were collected at 0 h, 15 min, 30 min, 1 h, 3 h, 6 h, 12 h, and 24 h after abiotic stress treatments.These materials were immediately frozen and stored in liquid nitrogen until total RNA extraction and real-time quantitative polymerase chain reaction (qRT-PCR).

RNA Isolation and Relative Expression Level Analysis of PheHsfs
For the tissue-specific expression analysis, the RPKM (the reads per kilobase of exon model per million mapped reads) values of PheHsfs were retrieved from transcriptome sequencing of developing flowers (four stages of flowering and leaves) and shoots of moso bamboo (winter shoots, six shoot heights, and culms) [1,46].RPKM values were used to analyze the relative expression levels of the PheHsf genes.For the flower development samples, four developmental stages (F1: the floral buds begun to form; F2: the floral organs gradually matured; F3: the flowers were in full blossom; and F4: the embryo formation stage) were defined based on the anatomical structure of floral organs by Gao et al. [46].For F1 and F2, the buds were collected, respectively.For F3 and F4, spikelets were collected, respectively.Leaves collected from non-flowering moso bamboo were defined as CK1.For the shoot development stage sample, four development stages of shoots were defined based on the continuing measurement of bamboo shoot height and anatomical changes by Li et al. [1].The eight samples according to the four developmental stages of moso bamboo shoots were as follows: S1 (winter shoots), S2-S5 (0.5 m, 1 m, 3 m, and 6 m, early growth period), S6-S7 (9 m and 12 m, late growth), and CK (clum after leaf expansion, mature period).For S1-S7, the shoot tips of different heights were collected, respectively.For CK2, each top internode was cut from the top to 1/2, then each top internode was divided into basal, middle and top.After that, the samples were cut from the tissue located in the top part of the three internodes above and then mixed.
For qRT-PCR, PheHsfs primers were designed using Primer 3.0 (http://primer3.ut.ee/).Primer sequences, amplicon Length, amplification efficiency, and correlation coefficients are listed in Table S1, and their specificity was verified using the online tool Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi) and the melting curves of PCR products.For every primer pair, a standard curve was constructed to calculate the gene-specific PCR efficiency from the 10-fold series dilution of the mix cDNA template.The R 2 (correlation coefficients) and slope values can be obtained from the standard curve.The following formula was used to calculate the corresponding PCR amplification efficiencies (E): E = (10 −1/slope − 1) × 100 [47].Tonoplast intrinsic protein 41 gene (TIP41) was used as an internal control [48].The qRT-PCR reactions were conducted using a SYBR Green I master mix (Roche, Mannheim, Germany).The qRT-PCR conditions were as follows: 45 cycles of 95 • C for 10 s, 60 • C for 10 s, and 72 • C for 20 s.Three replicates were performed for each gene.Gene expression was evaluated using the 2 −∆∆ Ct method [49].

Co-Expression Network and Protein Interactions of PheHsfs
The expression correlation of the PheHsfs was calculated by Pearson correlation coefficient (PCC; R-value) using gene expression RPKM values from the high-throughput transcriptome data in R. Expression correlation data were used for the correlation network, and co-expressed gene pairs were filtered with a PCC cut-off of 0.85 as previously described [50].Cytoscape version 3.4.0were used to analyze and visualize the network [51].
For the protein interaction networks, the homolog Hsf proteins in rice were constructed by STRING (http://stringdb.org/)using an option value >0.7.The homolog proteins of the determined interactive rice proteins were identified in moso bamboo by BLASTP analysis.

The Hsf Family Genes in Moso Bamboo
To identify Hsf genes in moso bamboo, a total of 41 candidate PheHsf genes were retrieved from the annotation in the Bamboo Genome Database.From these, 19 candidate PheHsf genes with incomplete HSF domains were considered as PheHsf-like genes, which were not selected for subsequent analysis in this study.Twenty-two putative moso bamboo PheHsf genes containing full HSF domains (PF00447) were confirmed by searching the Pfam and the Conserved Domain Databases.However, the CDS sequence for PheHsfA5, PheHsfB4a-1 and PheHsfB4c-1 contained 42 bases (from 268 to 309 base), 102 bases (233 to 336) and 123 base (256 to 377) inserts, respectively, when compared to the cloned cDNA sequence (Fasta S3).For this analysis, the CDS and amino acid sequences of these three genes are based on the cDNA sequences.The CDS and amino acid sequences are listed in supplementary files (Fasta S1 and S2).The identified 22 PheHsf genes were distributed among 22 scaffolds (Table 1).The sequences of 22 PheHsf genes were named according to the corresponding relationship among in rice, B. distachyon and moso bamboo.The number of amino acid (aa) sequences of PheHsf proteins ranged from 247 (PheHsfC1b-1) to 679 (PheHsfA1a), the predicted isoelectric points (pI) varied between 4.70 (PheHsfA6a) and 9.81 (PheHsfB4c-1), and the molecular weight (MW) ranged from 26.77 kDa (PheHsfC1b-1) to 75.05 kDa (PheHsfA1a) (Table 1).

Structure and Motif Analyses of PheHsf Genes
To better understand the gene structure diversity of PheHsfs, we compared the intron-exon arrangements and the conserved motifs (Figure S1).The number of introns in the Hsf genes of moso bamboo ranged from zero to three.Most of the PheHsfs (19/22) contained one to two introns (Figure S1b).Three introns were found in PheHsf4, whereas none were detected in PheHsfB4a-1 and PheHsfB4c-1.
Based on the known information on functional domains of Hsfs in some model plants [15,54], the sequence and positions of similar domains were identified in the PheHsfs by sequence alignment.As shown in Table 2, five conserved domains (DBD, HR-A/B, NLS, NES, and AHA) were identified.The DBD domain comprised of three α-helices (α1-3) and four β-sheets (β1-4) were found in all PheHsfs (Figure S2a).HR-A/B domain is critical for one Hsf interacting with other Hsfs to form trimer [15].All class A PheHsfs have a 21 amino-acid (aa) insertion between HR-A and HR-B regions; class C PheHsfs have seven aa insertions; all class B PheHsfs have no insertion (Figure S2b).NLS and NES domain function in the assembly of a nuclear import complex and the receptor-mediated export in complex with the NES receptor [15].The majority of the PheHsfs showed the presence of a NES and/or NLS domain.In addition, the activation domain AHA was found in all class A members but not in classes B and C (Table 2).

Structure and Motif Analyses of PheHsf Genes
To better understand the gene structure diversity of PheHsfs, we compared the intron-exon arrangements and the conserved motifs (Figure S1).The number of introns in the Hsf genes of moso bamboo ranged from zero to three.Most of the PheHsfs (19/22) contained one to two introns (Figure S1b).Three introns were found in PheHsf4, whereas none were detected in PheHsfB4a-1 and PheHsfB4c-1.
Based on the known information on functional domains of Hsfs in some model plants [15,54], the sequence and positions of similar domains were identified in the PheHsfs by sequence alignment.As shown in Table 2, five conserved domains (DBD, HR-A/B, NLS, NES, and AHA) were identified.The DBD domain comprised of three α-helices (α1-3) and four β-sheets (β1-4) were found in all PheHsfs (Figure S2a).HR-A/B domain is critical for one Hsf interacting with other Hsfs to form trimer [15].All class A PheHsfs have a 21 amino-acid (aa) insertion between HR-A and HR-B regions; class C PheHsfs have seven aa insertions; all class B PheHsfs have no insertion (Figure S2b).NLS and NES

Cis-Regulatory Element Analysis in Promoters of PheHsfs
To predict the biological function of PheHsfs, the 1500 bp upstream sequence from the translation start sites of PheHsf genes were analyzed using the PlantCARE database.The results show that the promoter of each PheHsf has several cis-regulatory elements such as phytohormone-(abscisic acid, jasmonic acid and gibberellic acid), abiotic stress-(low temperature, heat stress, drought, and fungal elicitor), and developmental process-related elements.Figure S3 shows that the ABA-responsive element (ABRE), the MeJA-responsive element (CGTCA-motif), and SA-responsive element (TCA-element) were found in the promoters of 17, 16, and 11 PheHsf genes, respectively.The promoters of 12 and 10 PheHsf genes contained the HSE and LTR, respectively.MYB-binding sites involved in drought inducibility (MBS), fungal elicitor-responsive elements (Box-W1) and defense-and stress-responsive elements (TC-rich) were found in 17, 15 and 8 PheHsf genes, respectively.Additionally, meristem expression (CAT-box), meristem-specific activation (CCGTCC-box), and endosperm expression (Skn-1_motif) motifs were found in the 13, 10 and 18 PheHsf genes, respectively.These findings indicate that PheHsfs might be associated with different transcriptional regulatory mechanisms for developmental, hormone and stress processes.

Expression Pattern of the PheHsf Genes in Shoot and Flower Development
Based on the RNA-Seq data of different flowering developmental stages [43] and the internodes of shoots at different heights [1], a heat map was constructed according to the RPKM of 22 PheHsfs (Figure 2).During four flowering developmental stages of moso bamboo, PheHsf genes could be classified into four groups (A, B, C, and D) according to their relative expression levels (Figure 2a).Most of the PheHsfs were highly expressed (RPKM > 10) in at least one stage, and only four PheHsf genes were expressed at low levels (RPKM < 2) in at least two stages during floral development (Table S2).Four members (PheHsfB4a-1, PheHsf15, PheHsf12, and PheHsfB4c-1) of group C had higher transcript accumulation during two earlier stages (F1 and F2) and were downregulated at two later stages (F3 and F4 stage).Group D consisted of 10 PheHsf genes (PheHsfA4a-1, PheHsfB2c, PheHsfA1a, PheHsfA4a-2, PheHsfC2b, PheHsfA7a, PheHsfB2a, PheHsfA2a-2, PheHsfA6b-1, and PheHsfB1), and their expression levels steadily increased at the F3 and F4 stages.Group B comprised of six PheHsfs (PheHsfA4d-1, PheHsfA4d-2, PheHsfA5, PheHsfA6b-2, PheHsfA6a, and PheHsfA2a-1), showing higher transcript accumulation at two later stages (F3 and F4) and in the leaves.Group A only had two genes, PheHsfC1b-1 and PheHsfC1b-2, with expression levels six times and three times higher in leaves than the four stages of floral development, respectively.During shoot growth, Most PheHsfs (19/22) Forests 2019, 10, 100 9 of 17 showed very low expression levels (RPKM < 8) at all seven stages, and only two genes (PheHsfA1a and PheHsfA6b-2) had higher expression levels (RPKM > 20) in at least one stage of bamboo shoot growth (Table S3).Based on the expression profiles, the PheHsfs were classified into four groups (Figure 2b).Among these, PheHsfA1a and PheHsfA2a-2 were clustered into the same group with continuous down-regulated expression from the S1 to S7 stage.However, PheHsfA6b-2 and PheHsfB2c were clustered in the same group, which showed twin peaks during the S2 and S5 stage, respectively.
PheHsfA4a-2, PheHsfC2b, PheHsfA7a, PheHsfB2a, PheHsfA2a-2, PheHsfA6b-1, and PheHsfB1), and their expression levels steadily increased at the F3 and F4 stages.Group B comprised of six PheHsfs (PheHsfA4d-1, PheHsfA4d-2, PheHsfA5, PheHsfA6b-2, PheHsfA6a, and PheHsfA2a-1), showing higher transcript accumulation at two later stages (F3 and F4) and in the leaves.Group A only had two genes, PheHsfC1b-1 and PheHsfC1b-2, with expression levels six times and three times higher in leaves than the four stages of floral development, respectively.During shoot growth, Most PheHsfs (19/22) showed very low expression levels (RPKM < 8) at all seven stages, and only two genes (PheHsfA1a and PheHsfA6b-2) had higher expression levels (RPKM > 20) in at least one stage of bamboo shoot growth (Table S3).Based on the expression profiles, the PheHsfs were classified into four groups (Figure 2b).Among these, PheHsfA1a and PheHsfA2a-2 were clustered into the same group with continuous down-regulated expression from the S1 to S7 stage.However, PheHsfA6b-2 and PheHsfB2c were clustered in the same group, which showed twin peaks during the S2 and S5 stage, respectively.S2 and S3.
For different periods of flower and shoot development in moso bamboo, the PheHsfA1a and PheHsfA6b-2 genes showed high transcript accumulation in all 11 stages, whereas four genes (PheHsfB4a-1, PheHsfA4d-1, PheHsfA7a, and PheHsfB4a-2) exhibited extremely low transcript accumulation at all stages.Moreover, the other 15 PheHsfs depicted significantly higher transcript levels during flower development than during shoot growth.These findings indicate that these genes have different regulatory roles in moso bamboo flower development and shoot growth.S2 and S3.
For different periods of flower and shoot development in moso bamboo, the PheHsfA1a and PheHsfA6b-2 genes showed high transcript accumulation in all 11 stages, whereas four genes (PheHsfB4a-1, PheHsfA4d-1, PheHsfA7a, and PheHsfB4a-2) exhibited extremely low transcript accumulation at all stages.Moreover, the other 15 PheHsfs depicted significantly higher transcript levels during flower development than during shoot growth.These findings indicate that these genes have different regulatory roles in moso bamboo flower development and shoot growth.

Expression Correlation and Interaction Networks
To further investigate the PheHsf proteins and how they interact with each other, a coexpression network was constructed using expression values of PheHsf genes during shoot and flower development.The connecting gene with PCC magnitude >0.85 was recognized as strongly coexpressed genes [36].The result showed that 18 PheHsf genes and 38 correlations in a co-expression network with PPC >0.85 cutoffs were obtained (Figure 4a).Genes with stronger correlation might play roles as interacting partners in similar biological pathways.Based on the results, PheHsfA4a-1 and PheHsfA4a-2 were recognized as hub genes as nodes with 10 and 9 connectivity in the whole network, respectively, which had more connectivity than other PheHsf genes.Both PheHsfA4a-1 and PheHsfA4a-2 had an up-regulated expression in response to different stress treatments (heat stress, cold, drought and salt).These results indicated that the two PheHsfA4s with a greater role in shoot and flower development also had an important role in stress responses.S4.

Expression Correlation and Interaction Networks
To further investigate the PheHsf proteins and how they interact with each other, a co-expression network was constructed using expression values of PheHsf genes during shoot and flower development.The connecting gene with PCC magnitude >0.85 was recognized as strongly co-expressed genes [36].The result showed that 18 PheHsf genes and 38 correlations in a co-expression network with PPC >0.85 cutoffs were obtained (Figure 4a).Genes with stronger correlation might play roles as interacting partners in similar biological pathways.Based on the results, PheHsfA4a-1 and PheHsfA4a-2 were recognized as hub genes as nodes with 10 and 9 connectivity in the whole network, respectively, which had more connectivity than other PheHsf genes.Both PheHsfA4a-1 and PheHsfA4a-2 had an up-regulated expression in response to different stress treatments (heat stress, cold, drought and salt).These results indicated that the two PheHsfA4s with a greater role in shoot and flower development also had an important role in stress responses.
To identify the two PheHsfA4a-associated proteins and protein complexes, prediction networks were built with STRING (http://www.string-db.org/)based on the interaction network of rice orthologous genes (Figure 4b).Because the rice orthologous gene of both PheHsfA4a-1 and PheHsfA4a-2 were OsHsfA4a, the identified moso bamboo proteins predicted to participate in the interaction network with PheHsfA4a-1 and PheHsfA4a-2 were the same.They interacted directly with 10 identified proteins, including five HSP70 proteins, one HSP90 protein, and four MAPK proteins.Because PheHsfA2a-2 (~22-fold higher than the control), PheHsfA4d-2 (~20-fold higher than the control) and PheHsfA4d-1 (~79-fold higher than the control) were strongly induced by heat, drought and NaCl stress, respectively.We also identified the two PheHsfA4d-and PheHsfA2a-2associated proteins and protein complexes, and PheHsfA4ds and PheHsfA2a-2 were also predicted to interact with HSP70 proteins, HSP90 protein, and MAPK proteins (Figure 4c,d).

Characterization of the Moso Bamboo Hsf Genes Family
Hsf genes play essential roles in plant development and in responding to various stress conditions [55,56].To explore the functions of PheHsf s in moso bamboo for the first time, the present study identified a total of 22 PheHsf genes according to the moso bamboo genome database and FLNC reads database [32,34,57].
We found that the moso bamboo has a similar number of Hsf s as rice, B. distachyon, maize, and A. thaliana (22)(23)(24)(25).This partially accounts for the support of Hsf conservation during evolution.Phylogenetic analysis of Hsfs in moso bamboo, O. sativa, B. distachyon, and A. thaliana indicated that PheHsfs have a higher degree of sequence similarity with OsHsfs and BdHsfs than AtHsfs, which coincides with the evolutionary relationships among the four species.All three Hsf classes (A, B, and C) were identified in the three monocots and one dicot, implying that the Hsf genes originated prior to the divergence of monocots and dicots.
In the investigation of conserved Hsf domains and intron-exon structures, all 22 of the PheHsfs contain the necessary (DBD and OD) and/or specific protein domains (NLS, NES, and AHA).The hydrophobic core of DBD domain ensures precise positioning and highly selective interaction with heat stress promoter elements [15].The OD of plant Hsfs confers distinct patterns of specificity for hetero oligomerization [15].These AHA motifs, which are located in the C-terminal, are characterized by aromatic (W, F, Y), large hydrophobic (L, I, V) and acidic (E, D) amino acid residues [15].Those domains might be essential for functional conservation [56].Twenty of twenty-two PheHsf genes have one intron in their DBD domain (Figure S1), which is an evolutionarily conserved intron [13].However, no intron was identified in PheHsfB4a-1 and PheHsfB4c-1 of subclass B4, which is different from subclass B4 of rice and B. distachyon [53].

Cis-Regulatory Element Analysis in the Promoters of PheHsfs
Previous studies have illustrated the key roles of Hsf s in developmental processes and stress tolerance through their regulation of target genes [58].The number and form of cis-elements in the promoter region might play an essential function in the regulation of gene expression in relation to metabolic pathways [59].The in silico survey of putative cis-elements using PlantCARE showed that 12 of the 22 PheHsfs promoters contained HSEs.This implies that PheHsfs might regulate themselves [55].Additionally, the promoter region of PheHsfA1a contained more types of development-related cis-elements (Figure S3), which coincides with its constitutive expression during shoot and flower development in moso bamboo (Figure 2).

PheHsfAs Involvement in Development Processes
The 22 PheHsfs exhibited diverse expression patterns during shoot and flower development of moso bamboo under normal conditions.PheHsfA1a was upregulated during shoot and flower development and was constitutively expressed in different tissues (Figure 2).These findings are similar to that in Arabidopsis [60] and Salix suchowensis [61].Under normal conditions, class A1 Hsfs of Arabidopsis are involved in housekeeping processes, and SsuHsf-A1a of Salix suchowensis are constitutively expressed in different tissues [60,61].In rice, nearly all classA members showed high expression levels in all tissues [45,62].In this study, nearly all the PheHsfAs (except PheHsfA7a) were found to show high transcription levels in the leaves, culms, and the four stages of flower development, similar to that in rice.Only two members of PheHsfAs (PheHsfA1a and PheHsfA6b-2) have high transcription levels in the seven stages of shoot development.This indicated that the function of PheHsfAs is conserved and/or specific in regulating flower development and shoot growth in moso bamboo.

PheHsfAs are Involved in Stress Responses
Under heat or other stress conditions, plant Hsfs regulate the transcription of target genes (Hsps and other stress-inducible genes) to enhance stress resistance.Recent genome-wide expression profile analyses showed that most of the Hsf genes are upregulated after heat, cold, drought, and salt stress [56,58,60].An increase in the number of Hsf genes has been shown to improve plant stress tolerance [54,63].In moso bamboo, Zhao et al. identified seven and two PheHsfs that were upregulated after 0.5 h and 8 h of high light stress (1200 µmol • m −2 • s −1 ), suggesting that these play vital roles both in response to short-term (0.5 h) and mid-term (8 h) high light [64].However, the regulatory roles of Hsfs in response to other abiotic stresses are unclear.Classes B and C do not harbor AHA motifs, which are essential for the activity of class A Hsfs [16].Therefore, class A PheHsf genes were selected to further study their response patterns under stress and hormone treatments.
Under heat treatment, 9 of 12 OsHsfAs and 9 of 13 BdHsfAs were upregulated in rice and B. distachyon, respectively [59,62].In our study, 8 of 12 PheHsfAs were found to be induced by high temperature.Of these, PheHsfA4d-1, PheHsfA2a-2 and PheHsfA6b responded more quickly to heat stress than others.PheHsfA2a-2 showed stronger induced under-heat stress conditions, which is similar to that of its homologous genes in rice, Arabidopsis, and B. distachyon [18,45,53].These findings suggest that PheHsfA2a-2 plays an important role in response to heat stress of moso bamboo.However, PheHsfA2a-1, which is the most similar to PheHsfA2a-2 and also belongs to A2a-type PheHsfs, was slightly upregulated 1 h after heat stress application.PheHsfA4d-2, which is the most similar to PheHsfA4d-1, was less upregulated by heat and salt stress compared to PheHsfA4d-1, but showed higher relative expression levels than PheHsfA4d-1 with cold and drought stresses.

Expression Correlation and Interaction Networks
In addition, the expression values of PheHsf genes in shoot and flower developmental stages were used to identify potential underlying co-expression networks using Cytoscape.Based on their degree of connectivity, the hub genes PheHsfA4a-1 and PheHsfA4a-2 were identified to play a regulatory role and correlated with other PheHsfs in the complex feedback network.According to the prediction of five PheHsf proteins (four PheHsfA4s and one PheHsfA2a-2) interaction networks, the interactive proteins might include MAPKs and heat shock proteins.MPK5 protein in rice acts as a positive regulator of drought, salt and cold tolerance; is involved in disease resistance and abiotic stress tolerance signaling pathways; and also negatively modulates pathogenesis-related (PR) gene expression and broad-spectrum disease resistance [66].The MPK1 protein in rice acts downstream of heterotrimeric G protein alpha subunit and small GTPase RAC1 and may regulate the expression of various genes involved in biotic and abiotic stress response [67].Base on the above information and our results, four PheHsfA4s and one PheHsfA2a-2 may play a very important role in shoot and flower development and stress response.

Conclusions
In this study, 22 PheHsf genes in moso bamboo were identified for the first time.These genes could be classified into three classes (A, B and C) according to the comparison of the phylogenetic relationships with O. sativa, B. distachyon and A. thaliana Hsf genes.Expression analyses revealed that two hub genes, PheHsfA4a-1 and PheHsfA4a-2, might act as a potential "node" for crosstalk between developmental processes and abiotic stress responses.Furthermore, PheHsfA2a-2, PheHsfA4d-1 and PheHsfA4d-2 might also act as essential parts in response to stress.These results provide insights into the responses of PheHsfAs to abiotic stresses treatments, although their underlying molecular mechanism requires further study.
Author Contributions: X.L. and J.G. designed the experiments; X.L. and J.L. (Juan Li) performed the tissue and organ collection; L.X. and X.L. performed the experiments; Z.C. formal analysis; J.L. (Jun Liu) and S.M. assisted in equipment maintenance and sample collection; L.X. writing-original draft preparation; X.L. and D.H. writing-review and editing; J.G. review and funding acquisition.L.X. and X.L. contributed equally in this work and should be considered first co-authors.

Figure 1 .
Figure 1.Phylogenetic relationship of PheHsf, AtHsf and OsHsf proteins.Neighbor-joining method and MEGA 7.0 software were used for phylogenetic analysis of Hsf proteins from Phyllostachys edulis (22 PheHsfs), Arabidopsis thaliana (23 AtHsfs), Oryza sativa (25 OsHsfs), and Brachypodium distachyon (24 BdHsfs).The names of subclass are shown outside of the circle.Branch lines of subclass are colored, indicating different Hsf subclasses.A MEME motif search revealed a total of 19 motifs in the PheHsfs (FigureS1c).Three motifs (1, 2, and 4) constituting the DBD domain were identified.Motif 3 and 5 for the OD domain were detected in class A and C, and motif 6 for the OD domain was observed in class B. Motif 6 and motif 11, motif 8, and motif 17 (NSL) were observed in class A, class B and C, respectively.In general, the structure of the PheHsf proteins is conserved in moso bamboo.Furthermore, motif 7 represented the AHA domain close to the PheHsfA C-terminus (FigureS1cand Table2).

Figure 2 .
Figure 2. Expression pattern of PheHsfs in developmental flowers and shoots of moso bamboo.(a) The expression profile of PheHsfs in different stages of flowering.F1-F4: The four stages of developmental flowers (the floral buds begun to form, the floral organs gradually matured but did not undergo flowering, the flowers were in full blossom, and the embryo formation stage); CK1: leaves.(b) The expression profile of PheHsfs in different stages of shoots.S1: winter shoot; S2-S7: different heights of shoots (0.5 m, 1 m, 3 m, 6 m, 9 m and 12 m), and CK2 (culms after leaf expansion).All the samples had one repeat.The heatmap was pictured using R based on the RPKM of Phehsfs in each sample.Details of the RPKM are shown in TablesS2 and S3.

Figure 2 .
Figure 2. Expression pattern of PheHsfs in developmental flowers and shoots of moso bamboo.(a) The expression profile of PheHsfs in different stages of flowering.F1-F4: The four stages of developmental flowers (the floral buds begun to form, the floral organs gradually matured but did not undergo flowering, the flowers were in full blossom, and the embryo formation stage); CK1: leaves.(b) The expression profile of PheHsfs in different stages of shoots.S1: winter shoot; S2-S7: different heights of shoots (0.5 m, 1 m, 3 m, 6 m, 9 m and 12 m), and CK2 (culms after leaf expansion).All the samples had one repeat.The heatmap was pictured using R based on the RPKM of Phehsfs in each sample.Details of the RPKM are shown in TablesS2 and S3.

Figure 3 .
Figure 3. Expression analysis of PheHsfAs under different abiotic stress treatments in moso bamboo.(a) Heat map representation for the expression patterns of PheHsfAs after 15 min, 30 min, 1 h, 3 h, 6 h, 12 h, and 24 h of heat, cold, drought, and salt stresses: expression levels under stress vs. control; the expression results were obtained by qRT-PCR.The different colors correspond to log2 transformed value.(b) Venn diagram showing the number of overlapping PheHsfAs that are upregulated > two-fold under abiotic stress: heat, cold, salt, and drought.Details of the expression data are shown in TableS4.

Figure 3 .
Figure 3. Expression analysis of PheHsfAs under different abiotic stress treatments in moso bamboo.(a) Heat map representation for the expression patterns of PheHsfAs after 15 min, 30 min, 1 h, 3 h, 6 h, 12 h, and 24 h of heat, cold, drought, and salt stresses: expression levels under stress vs. control; the expression results were obtained by qRT-PCR.The different colors correspond to log2 transformed value.(b) Venn diagram showing the number of overlapping PheHsfAs that are up-regulated > two-fold under abiotic stress: heat, cold, salt, and drought.Details of the expression data are shown in TableS4.

1 cFigure 4 .
Figure 4. Co-expression network and interaction network of selected PheHsf genes in moso bamboo.(a) The model was built based on RPKM of RNA-seq of shoot and flower development.(b-d) Interaction network of PheHsfA4as, PheHsfA4ds and PheHsfA2a-2 in moso bamboo, respectively.Colored balls (protein nodes) in the network were used as a visual aid to indicate different input proteins and predicted interactors.Protein nodes which are enlarged indicate the availability of 3D protein structure information.Gray lines connect proteins which are associated by recurring textmining evidence.Line thickness indicates the strength of data support.

Table 1 .
Overview of PheHsf genes in moso bamboo.