Genome-Wide Identification and Expression Analysis of Heat Shock Transcription Factor Family in Chenopodium quinoa Willd

Heat shock transcription factor (Hsf) is one of the conserved gene families in plants, playing a crucial role in growth and development, as well as in response to diverse stresses. Although it has been systematically studied in many species, little is known about the Hsf gene family in Chenopodium quinoa, especially those involved in the regulatory network of stress processes. In this study, we identified 23 Hsf genes in quinoa (CqHsfs) through a genome-wide search method based on the latest available genome information. Phylogenetic analysis classified them into three groups, and group A was further divided into nine subgroups, which was supported by conserved domain organizations. Gene structure and multiple sequence alignment analysis revealed that all of the CqHsfs possessed a similar structure organization and were highly conserved in BDB domain. Interaction network analysis identified 13 CqHsfs involved in the network pathway to regulate diverse biological processes. Expression profiles of these CqHsfs were further investigated using the RNA-seq data, and tissue-specific and stress-responsive candidates were identified. Finally, four heat-responsive CqHsfs were selected to validate their expression level through semi-quantitative RT-PCR analysis. This study reported the organization, structure, and expression profiles of the Hsf gene family in quinoa, which will contributes to further functional analysis, and helps to better understand the roles and regulatory mechanism of heat shock factors playing in quinoa and beyond.


Introduction
Salinity, drought, heat, and cold as well as chemical toxicity are the main environmental stresses that damage plant growth and development to cause serious yield reduction of crop productivity and economic loss annually [1]. Among them, heat stress has gradually become the most detrimental factors with the advent of global warming. Generally, high temperatures lead to cellular damage and cell program death, resulting in the impairment of plant growth and development throughout the whole ontogeny, such as seed germination, photosynthesis, respiration, and membrane stability [2]. In addition, heat stress induces the excess reactive oxygen species (ROS), which leads to oxidative stress. In the long-term evolutional process, plants have implemented various mechanisms to cope with heat stress, of which heat shock transcription factors (Hsfs) acted as the crucial factor to regulate this process. Extensive studies has revealed that Hsfs were the terminal of the heat signal transduction

Identification of Hsf Genes in the Quinoa Genome
The potential members of the quinoa Hsf gene family were identified following the method described by Wang et al. [25] with some modifications. Firstly, all the available quinoa protein sequences were obtained from the Ensembl Plants database as a local protein database [26]. Then, BLASTP (Basic Local Alignment Search Tool Proteins) search was performed using the Hsf genes of arabidopsis as queries with the e-value of 1 × 10 −5 and identity of 50% as the threshold. Furthermore, the available hidden Markov model (HMM) profile of the Hsf gene family was downloaded from the PFAM database [27], and then an HMM search was performed using the hmmsearch tools embedded in HMMER3.0 [28]. Moreover, the available genome and protein sequences of quinoa, which were annotated as Hsf genes, were also downloaded from the NCBI (National Center for Biotechnology Information) database (https://www.ncbi.nlm.nih.gov). After manual redundant removal, the remaining sequences were submitted to NCBI-CDD (Conserved domain database) (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and the PFAM database to identify the presence and integrity of the conserved protein domain to obtain the candidate Hsf genes. The Mw (molecular weight), amino acids length, and theoretical pI (isoelectric point) of these identified Hsf genes were predicted by ExPASy online software.

Multiple Alignment, Phylogenetic, and Conserved Motif Analysis of CqHsfs
The Hsf protein sequences of four dicot species, including Arabidopsis (Arabidopsis thaliana L.), Poplar (Populus trichocarpa (Torr. & Gray)), soybean (Glycine max L.), and rape (Brassica rapa L.) together with four monocot species, including rice (Oryza sativa), Brachypodium (Brachypodium distachyon L.), barley (Hordeum vulgare L.), as well as maize (Zea mays L.) were downloaded from Heatsterdatabase (http://www.cibiv.at/services/hsf/db/search) [16]. Multiple sequence alignments were performed using ClustalW [29], then the phylogenetic tree was constructed using MEGA7.0 software [30] with the neighbor-joining (NJ) method, and the bootstrap was set to 1000 replications. The conserved protein domains of CqHsfs were predicted by the MEME online program [31]. The gene structure information was obtained from genome annotation file and displayed by the Gene Structure Display Server program [32].

Interaction Network Analysis of CqHsfs
The interaction network of CqHsfs was constructed based on the homologous genes between quinoa and arabidopsis by submitting the CpHsf genes to INPARANOID8 (http://inparanoid.sbc. su.se/cgi-bin/index.cgi). Then, AraNetV2 (http://www.inetbio.org/aranet) was used to obtain the homologous pairs with arabidopsis (AtHsf). Finally, the interaction network of CqHsfs was visualized by using Cytoscape_v3.2.1 software (http://www.cytoscape.org/) with the edge information.

Expression Profiles Analysis of CqHsf Genes
The available RNA-seq data of quinoa in different organs and under different stresses were obtained from the NCBI Sequence Read Archive (SRA) database (https://www.ncbi.nlm.nih.gov/sra). Detail information of the dataset we used are listed in Table S1. HISAT2 and StringTie software were used to analyze gene expression level [33]. The FPKM value (fragments per kilobase of transcript per million fragments mapped) of each gene was calculated, and the log10-transformed (FPKM + 1) values of quinoa Hsf genes were used to generate heat maps by using the Pheatmap package in R v3.0.0 software (Lucent Technologies, Murray Hill, NJ, USA).

Plant Growth, Heat Stress Treatment, and RT-PCR Analysis
Seeds of quinoa cultivar "Real Blanca" were grown in a growth chamber at controlled conditions (20 ± 1 • C, 16-h light/8-h dark cycle). Two week old seedlings were used to perform heat stress treatment. Heat stress was performed by putting the quinoa seedlings into the growth chamber with a temperature of 37 • C. Then, the leaves were collected for RNA extraction at 0, 6, 12, and 24 h after treatment. Total RNA was isolated using RNAiso Reagent (TaKaRa, Dalian, China) according to the manufacturer's instructions. Then, 4 heat-responsive CqHsfs were randomly selected to validate the expression patterns based on semi-quantitative RT-PCR analysis follow the Chen et al.'s [34] method. The used primers are listed in Table S2 and actin (LOC110724665) was used as the inner reference. Three biological replicates for each sample were performed.

Genome-Wide Identification of Hsf Genes in Quinoa
To identify the potential members of Hsf genes in quinoa, both the local BLASTP and HMM search were used. After manual correction, a total of 23 unique genes with conserved HSF domains were identified as the Hsf genes in quinoa. It reported that Hsf family was one of the smallest families in angiosperms with 21 members in arabidopsis, 25 in rice, 19 in castor bean, and 27 in poplar [16]. The size of Hsf genes in quinoa is similar with that of arabidopsis and rice, but significantly lower that of soybean (52) in which recent whole genome duplication (WGD) events have occurred due to palaeopolyploid [16], suggesting there was no recent WGD or gene tandem and segmental duplication occurred in quinoa genome. Furthermore, the identified quinoa Hsf genes were named as CqHsf1 to CqHsf23 according to chromosome localization. The amino acids length, Mw, and theoretical pI of these Hsf genes were also calculated. As shown in Table 1, the protein length of putative CqHsf genes ranged from 251 to 546 amino acids, with molecular weight (Mw) ranged from 28.9 to 60.7 kDa and theoretical pI ranged from 4.74 to 8.49, respectively, showing that there were significant differences among CqHsfs in physical and chemical properties (Table 1).

Multiple Alignment, Phylogenetic and Conserved Domains Analysis of CqHsfs
In order to understand the grouping and phylogenetic relationships of Hsf genes in quinoa, all 23 CqHsf protein sequences, together with 239 Hsfs from eight other species were used to construct the phylogenetic tree. As shown in Figure 1, the Hsf genes from different species were mainly clustered into three groups, namely A, B, and C groups. Then, the A group could be further divided into nine subgroups, which is consistent with previous reports [16,20]. Furthermore, Hsf genes from different species were unevenly clustered into these groups, indicating that the parallel evolutionary event of Hsf genes were occurred among these species, and homologous genes of model plants could provide some clues for revealing the biological function of counterpart Hsfs in quinoa.

Multiple Alignment, Phylogenetic and Conserved Domains Analysis of CqHsfs
In order to understand the grouping and phylogenetic relationships of Hsf genes in quinoa, all 23 CqHsf protein sequences, together with 239 Hsfs from eight other species were used to construct the phylogenetic tree. As shown in Figure 1, the Hsf genes from different species were mainly clustered into three groups, namely A, B, and C groups. Then, the A group could be further divided into nine subgroups, which is consistent with previous reports [16,20]. Furthermore, Hsf genes from different species were unevenly clustered into these groups, indicating that the parallel evolutionary event of Hsf genes were occurred among these species, and homologous genes of model plants could provide some clues for revealing the biological function of counterpart Hsfs in quinoa. Then, the conserved protein motifs in these quinoa Hsf genes were predicted using the MEME online tool. A total of 15 motifs were identified ( Figure 2). All of the quinoa Hsf genes contained motif1, motif2, and motif3. In addition, some motifs were only found in a certain group of CqHsfs, such as motif5 and motif7 were found in group B specifically, and motif6 only existed in the group A2 subfamily. It showed that the same cluster of the CqHsfs by phylogenetic analysis shared the similar conserved domain composition, suggesting that the protein structure is remarkably conserved within the subgroup of CqHsfs, and the specific motifs might play a role in determining the specific expression or functions. Gene structure of CqHsfs was further displayed using the GSDS online tool. Result shows that all of the CqHsf genes have a similar gene structure with only two exons and one intron ( Figure 2B).
To better understand the structural features and potential functions of CqHsf gene family, the conserved functional domains of them were further investigated. A total of seven different Then, the conserved protein motifs in these quinoa Hsf genes were predicted using the MEME online tool. A total of 15 motifs were identified ( Figure 2). All of the quinoa Hsf genes contained motif1, motif2, and motif3. In addition, some motifs were only found in a certain group of CqHsfs, such as motif5 and motif7 were found in group B specifically, and motif6 only existed in the group A2 subfamily. It showed that the same cluster of the CqHsfs by phylogenetic analysis shared the similar conserved domain composition, suggesting that the protein structure is remarkably conserved within the subgroup of CqHsfs, and the specific motifs might play a role in determining the specific expression or functions. Gene structure of CqHsfs was further displayed using the GSDS online tool. Result shows that all of the CqHsf genes have a similar gene structure with only two exons and one intron ( Figure 2B). B, and C) according to the peculiarities of the flexible linkers, HR-A/B regions, and RD conserved domains [16,20]. Among the 23 identified CqHsfs, twelve CqHsf genes belong to class A, nine belong to class B, and only two belong to class C. Furthermore, group A was divided into nine subfamilies according to the classification criteria as described by Scharf et al. [16], which were consistent with the phylogenetic analysis.

Interaction Network of CqHsfs
The interaction network of CqHsfs was predicted using the orthologous genes in arabidopsis. A total of 13 CqHsf genes were found to have orthologous with AtHsf genes (Figure 3). The 13 CqHsf genes belong to three groups, of which eight belong to group A, four belong to group B, and only one belong to group C. This result shows that the 13 CqHsf genes have interaction with a total of 151 other functional genes (Table S4). Among them, CqHsf 13, CqHsf15, CqHsf19, and CqHsf22 interacted To better understand the structural features and potential functions of CqHsf gene family, the conserved functional domains of them were further investigated. A total of seven different conserved functional domains were identified (Table S3), including DBD, HR-A/B, RD, NLS, AHA, RD, and NES. Among them, DBD was the most abundant domain, which is located in the N-terminus of all the CqHsf proteins, and consisted of about 100 amino acids. The DBD functional domain possesses a DNA-binding function and provides a specific recognition of the heat stress promoter elements (HSE). Multiple alignment of the DBD region showed that DBD was highly conserved in all of the CqHsfs (Figure 3). The oligomerization domain (HR-A/B region) was also found in all CqHsf proteins. Interestingly, RD conserved functional domain was specifically found in group B CqHsfs. In plants, the Hsf gene family could be divided into three categories (Classes A, B, and C) according to the peculiarities of the flexible linkers, HR-A/B regions, and RD conserved domains [16,20]. Among the 23 identified CqHsfs, twelve CqHsf genes belong to class A, nine belong to class B, and only two belong to class C. Furthermore, group A was divided into nine subfamilies according to the classification criteria as described by Scharf et al. [16], which were consistent with the phylogenetic analysis.

Interaction Network of CqHsfs
The interaction network of CqHsfs was predicted using the orthologous genes in arabidopsis. A total of 13 CqHsf genes were found to have orthologous with AtHsf genes (Figure 3). The 13 CqHsf genes belong to three groups, of which eight belong to group A, four belong to group B, and only one belong to group C. This result shows that the 13 CqHsf genes have interaction with a total of 151 other functional genes (Table S4). Among them, CqHsf 13, CqHsf15, CqHsf19, and CqHsf22 interacted with stress-related and growth-related transcription factor genes, such as ERF, WRKY, and MYB genes. The interaction network analysis of CqHsfs found that they play an important role in quinoa growth and development as well as stress response, which also contributed to better understanding the regulatory pathways which CqHsfs involved in controlling these biological processes. with stress-related and growth-related transcription factor genes, such as ERF, WRKY, and MYB genes. The interaction network analysis of CqHsfs found that they play an important role in quinoa growth and development as well as stress response, which also contributed to better understanding the regulatory pathways which CqHsfs involved in controlling these biological processes.

Expression Patterns of CqHsfs in Different Tissues or Organs
Gene expression patterns are usually closely related to their functions [25]. To gain insight into the temporal and spatial expression patterns of CqHsf genes, RNA-seq data of different organs, including seedling, stem, leaf, inflorescence, and dry seeds were downloaded from the NCBI Sequence Read Archive (SRA) database, and then were used to analyze the expression profiles in tissues. Results found that most of CqHsf genes were expressed in the five tissues ( Figure 4). Among them, CqHsf10 and CqHsf17 in group B were found to be highly expressed in seeding, stem, leaf, and inflorescence, while showing relatively low expression in dry seeds. CqHsf13 and CqHsf14 belonging to the A4 subgroup had relatively high expressions in all five organs, which was also reported in other plant species [35]. For example, the AtHsfA1 in arabidopsis, and MdHsfA1a, MdHsfA1d, MdHsfB1a in apple as well as HsfA1 and HsfA8 genes in wheat were found to be constitutive expression in the tested organs [36][37][38]. The high expression of CqHsf genes in most tissues indicates that these heat shock transcription factors might be widely involved in the growth and development processes. At the same time, CqGsf2, CqHsf3, CqHsf4, CqHsf9, and CqHsf19 showed relatively low expressions in almost all of the five tested tissues, indicating they might not be involved in the quinoa organ development. Interestingly, CqHsf19 specifically expressed in dry seeds, indicating that CqHsf19 may have a specific function in seed development and maturation.

Expression Patterns of CqHsfs in Different Tissues or Organs
Gene expression patterns are usually closely related to their functions [25]. To gain insight into the temporal and spatial expression patterns of CqHsf genes, RNA-seq data of different organs, including seedling, stem, leaf, inflorescence, and dry seeds were downloaded from the NCBI Sequence Read Archive (SRA) database, and then were used to analyze the expression profiles in tissues. Results found that most of CqHsf genes were expressed in the five tissues ( Figure 4). Among them, CqHsf10 and CqHsf17 in group B were found to be highly expressed in seeding, stem, leaf, and inflorescence, while showing relatively low expression in dry seeds. CqHsf13 and CqHsf14 belonging to the A4 subgroup had relatively high expressions in all five organs, which was also reported in other plant species [35]. For example, the AtHsfA1 in arabidopsis, and MdHsfA1a, MdHsfA1d, MdHsfB1a in apple as well as HsfA1 and HsfA8 genes in wheat were found to be constitutive expression in the tested organs [36][37][38]. The high expression of CqHsf genes in most tissues indicates that these heat shock transcription factors might be widely involved in the growth and development processes. At the same time, CqGsf2, CqHsf3, CqHsf4, CqHsf9, and CqHsf19 showed relatively low expressions in almost all of the five tested tissues, indicating they might not be involved in the quinoa organ development. Interestingly, CqHsf19 specifically expressed in dry seeds, indicating that CqHsf19 may have a specific function in seed development and maturation.

Expression Patterns of CqHsf Genes under Abiotic Stress Conditions
It is well known that the expression levels of Hsf genes are significantly affected by heat and other stresses in plants [16,38]. In this study, we also analyze the expression profiles of CqHsf genes under various stress conditions to identify the stress-responsive candidates. Under heat stress, 22 out of 23 CqHsfs were found to show differential expression ( Figure 5A), of which CqGsf2, CqHsf3, CqHsf4, CqGsf9, and CqHsf10 showed upregulated expression, while CqGsf14, CqHsf16, and CqHsf23 had downregulated expression both in roots and leaves. The specific expression patterns of them were also found under heat stress. CqHsf5 was found to be upregulated expression in leaves, while downregulated in roots, and CqHsf8 was upregulated 5.4 times compared to control in roots, while there was no significant variation in leaves. These heat-induced CqHsfs provided the elite candidates for further functional study. Furthermore, the expression patterns of these CqHsfs under drought and salt stresses were also investigated ( Figure 5B). Overall, the average expression level of these CqHsfs under drought and salt were significantly lower than that of heat stress, although differential expressions were also identified. CqHsf5 and CqHsf12 in group C were significantly upregulated under drought stress treatment, indicating that heat shock transcription factors in group C played a key role in drought stress. CqHsf8 and CqHsf17 were not only related to drought resistance, but also associated with salt tolerance. Interestingly, the expression of CqHsf7 was upregulated in leaf and leaf with epidermal bladder cells with 100 mM sodium chloride treatment and under drought stress, but was downregulated in epidermal bladder cells of 100-mM sodium chloride treated plants, indicating that CqHsf7 may have a special function to increase the adaptability of plants in response to abiotic stresses. In addition, CqHsf18 and CqHsf21, had similar expression profiles in different abiotic stresses, indicating that they might have analogous physiological functions. However, there

Expression Patterns of CqHsf Genes under Abiotic Stress Conditions
It is well known that the expression levels of Hsf genes are significantly affected by heat and other stresses in plants [16,38]. In this study, we also analyze the expression profiles of CqHsf genes under various stress conditions to identify the stress-responsive candidates. Under heat stress, 22 out of 23 CqHsfs were found to show differential expression ( Figure 5A), of which CqGsf2, CqHsf3, CqHsf4, CqGsf9, and CqHsf10 showed upregulated expression, while CqGsf14, CqHsf16, and CqHsf23 had downregulated expression both in roots and leaves. The specific expression patterns of them were also found under heat stress. CqHsf5 was found to be upregulated expression in leaves, while downregulated in roots, and CqHsf8 was upregulated 5.4 times compared to control in roots, while there was no significant variation in leaves. These heat-induced CqHsfs provided the elite candidates for further functional study. Furthermore, the expression patterns of these CqHsfs under drought and salt stresses were also investigated ( Figure 5B). Overall, the average expression level of these CqHsfs under drought and salt were significantly lower than that of heat stress, although differential expressions were also identified. CqHsf5 and CqHsf12 in group C were significantly upregulated under drought stress treatment, indicating that heat shock transcription factors in group C played a key role in drought stress. CqHsf8 and CqHsf17 were not only related to drought resistance, but also associated with salt tolerance. Interestingly, the expression of CqHsf7 was upregulated in leaf and leaf with epidermal bladder cells with 100 mM sodium chloride treatment and under drought stress, but was downregulated in epidermal bladder cells of 100-mM sodium chloride treated plants, indicating that CqHsf7 may have a special function to increase the adaptability of plants in response to abiotic stresses. In addition, CqHsf18 and CqHsf21, had similar expression profiles in different abiotic stresses, indicating that they might have analogous physiological functions. However, there are also some CqHsf genes, such as CqHsf4, CqHsf6, and CqHsf19, with little or no altered expression levels under various abiotic stress. The different expression patterns of CqHsf genes indicates that the function divergent has occurred in these CqHsfs. The identified stress-specific members provided excellent candidates for further functional analysis to reveal the function and molecular mechanism of CqHsfs playing in stress response. are also some CqHsf genes, such as CqHsf4, CqHsf6, and CqHsf19, with little or no altered expression levels under various abiotic stress. The different expression patterns of CqHsf genes indicates that the function divergent has occurred in these CqHsfs. The identified stress-specific members provided excellent candidates for further functional analysis to reveal the function and molecular mechanism of CqHsfs playing in stress response.

Validation of the Expression of CqHsfs by Semi-Quantitative RT-PCR Analysis
To verify the expression of these identified CqHsfs, four genes showing upregulation in leaves by RNA-seq analysis were randomly selected to detect their expression levels in four time points under 37 °C treatments through semi-quantitative RT-PCR ( Figure 6). Results showed all of the four genes were expressed, although different expression patterns were also found. At 0 h time point, all of them did not show any expression. Then, they increased rapidly their expression with heat treatment to the high level, and then gradually reduced the expression, suggesting they were heat-induced genes. CqHsf3 and CqHsf9 showed higher expression level at 6 h, while CqHsf4 and CqHsf10 had higher expression at 12 h. On the whole, the RT-PCR results validated the expression trends of these CqHsfs, which provided an important gene resource for future functional study of Hsfs in quinoa and beyond.

Validation of the Expression of CqHsfs by Semi-Quantitative RT-PCR Analysis
To verify the expression of these identified CqHsfs, four genes showing upregulation in leaves by RNA-seq analysis were randomly selected to detect their expression levels in four time points under 37 • C treatments through semi-quantitative RT-PCR ( Figure 6). Results showed all of the four genes were expressed, although different expression patterns were also found. At 0 h time point, all of them did not show any expression. Then, they increased rapidly their expression with heat treatment to the high level, and then gradually reduced the expression, suggesting they were heat-induced genes. CqHsf3 and CqHsf9 showed higher expression level at 6 h, while CqHsf4 and CqHsf10 had higher expression at 12 h. On the whole, the RT-PCR results validated the expression trends of these CqHsfs, which provided an important gene resource for future functional study of Hsfs in quinoa and beyond.

Conclusions
By performing a genome-wide search, we identified 23 quinoa Hsf genes in this study. These CqHsfs could be classified into three groups and the same group shared a similar domain organization and gene structure. Interaction network analysis found that 13 CqHsfs participated in regulating 153 other quinoa functional genes to form the complex network pathway about diverse biological processes. Expression profiles of these CqHsfs were further investigated using the RNA-seq data, and the tissue-specific and stress-responsive candidates were identified. Finally, four heat-responsive CqHsfs were selected to validate their expression level through semi RT-PCR analysis. Our study for the first time reported the genomic organization, evolution, and expression of the quinoa heat shock transcription factor family, which not only provides candidates for future functional analysis, but also lays the foundation for revealing the molecular mechanism of Hsfs regulating quinoa development and stress response.
Supplementary Materials: The following are available online at www.mdpi.com/link, Table S1: The accession numbers and samples information of RNA-seq data used in this study, Table S2: The primers were used for semi-quantitative RT-PCR analysis, Table S3: Conserved functional domains of CqHsfs.
Author Contributions: G.T. and H.Z. (Haoshuang Zhan) performed all analysis and prepared the paper. G.X. contributed to data collection and gene family identification. C.X. and H.Z. (Hong Zhang) contributed to data analysis. X.N. conceived the study and revised the manuscript. W.J. provided experimental coordination and reviewed the manuscript. All authors read and approved the final manuscript.

Conclusions
By performing a genome-wide search, we identified 23 quinoa Hsf genes in this study. These CqHsfs could be classified into three groups and the same group shared a similar domain organization and gene structure. Interaction network analysis found that 13 CqHsfs participated in regulating 153 other quinoa functional genes to form the complex network pathway about diverse biological processes. Expression profiles of these CqHsfs were further investigated using the RNA-seq data, and the tissue-specific and stress-responsive candidates were identified. Finally, four heat-responsive CqHsfs were selected to validate their expression level through semi RT-PCR analysis. Our study for the first time reported the genomic organization, evolution, and expression of the quinoa heat shock transcription factor family, which not only provides candidates for future functional analysis, but also lays the foundation for revealing the molecular mechanism of Hsfs regulating quinoa development and stress response.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4395/8/7/103/ s1, Table S1: The accession numbers and samples information of RNA-seq data used in this study, Table S2: The primers were used for semi-quantitative RT-PCR analysis, Table S3: Conserved functional domains of CqHsfs.
Author Contributions: G.T. and H.Z. (Haoshuang Zhan) performed all analysis and prepared the paper. G.X. contributed to data collection and gene family identification. X.C. and H.Z. (Hong Zhang) contributed to data analysis. X.N. conceived the study and revised the manuscript. W.J. provided experimental coordination and reviewed the manuscript. All authors read and approved the final manuscript.