HSP70 Gene Family in Brassica rapa: Genome-Wide Identification, Characterization, and Expression Patterns in Response to Heat and Cold Stress

Heat shock proteins protect plants from abiotic stress, such as salt, drought, heat, and cold stress. HSP70 is one of the major members of the heat shock protein family. To explore the mechanism of HSP70 in Brassica rapa, we identified 28 putative HSP70 gene family members using state-of-the-art bioinformatics-based tools and methods. Based on chromosomal mapping, HSP70 genes were the most differentially distributed on chromosome A03 and the least distributed on chromosome A05. Ka/Ks analysis revealed that B. rapa evolution was subjected to intense purifying selection of the HSP70 gene family. RNA-sequencing data and expression profiling showed that heat and cold stress induced HSP70 genes. The qRT-PCR results verified that the HSP70 genes in Chinese cabbage (Brassica rapa ssp. pekinensis) are stress-inducible under both cold and heat stress. The upregulated expression pattern of these genes indicated the potential of HSP70 to mitigate environmental stress. These findings further explain the molecular mechanism underlying the responses of HSP70 to heat and cold stress.


Introduction
The heat shock transcription factor family strengthens plants under biotic and abiotic stress conditions; as a result, plants proceed with normal growth and development [1]. Temperature significantly affects plant growth; all plants require an optimum temperature for survival. Thus, many studies have been conducted on plants to explore the response to temperature stress [2,3]. (Heat shock treatment has been shown to stimulate a genetic network of proteins as a defense system in plants, and these proteins are named heat shock proteins (HSPs) [4]. HSPs are divided into subgroups on the basis of their molecular weight (60, 70, 90, and 100 KDa), named HSP60, HSP70, HSP90, and HSP100 gene families, respectively [5]. Of these heat shock families, HSP70 has been identified as the most conserved throughout evolution and play vital role in the maintenance of normal plant development under heat stress that can be induced by exposure to sublethal temperatures in eukaryotes [6][7][8]. HSP70 act as chaperons to protect the plant under temperature stress condition by covering the heat and cold sensitive molecules. It helps to reform their folding structure and repairing the deformed structures. All the HSP70 proteins consists of two functional domains named as nucleotide binding domain (NBD) and substrate binding domain (SBD) [9].
The HSP70 family, which encodes HSP70 proteins, has been identified and studied in many plants, such as Arabidopsis (17 genes), rice (26 genes), spinach (12 genes), and soybean (61 genes) [8]. In Arabidopsis, HSP70 deficient plants showed the stunted plant growth and abnormal leaf phenotype [10], double HSP70 knockout leads to the defective male and female gametes [11]. Furthermore, the cytosolic and nuclear HSP70 proteins showed significant role in Arabidopsis plant development under abiotic stress conditions [12]. To date, little information of HSP70 genes has been elucidated in Brassica crops. In the current study, the HSP70 gene family will be explored in a Brassica rapa family member. This vegetable-based family includes cabbage, Chinese cabbage, mustard, and turnip, which are sources of oil, amino acids, and proteins for humans and animals. These vegetables are highly susceptible to temperature stress [13] and are adversely affected by heat, cool, drought, and biotic stress during their cultivation, which causes germination delays, growth arrest, frost damage, and poor head yield [14].
Brassica rapa has gained remarkable economic and nutritional value due to its heading trait, which makes it highly susceptible to stress. Under temperature stress, the leaves become pale and weakened due to chloroplast damage, and the plant becomes unable to properly regulate photosynthesis [15]. Under such conditions, the HSP70 gene family plays an important role in protecting the plant and reestablishing the cellular and biological functioning of plants. Genome-wide identification and characterization of the HSP70 gene family have not been studied in Brassica rapa to date. We identified and explored the members of the HSP70 gene family, constructed their evolutionary tree and synteny correlations, and predicted protein structures with available bioinformatics software. Moreover, we analyzed the expression of these proteins by qRT-PCR to determine their responses to heat shock stress in Chinese cabbage (Brassica rapa subsp. Pekinensis). These findings will open up new gates for understanding the structure, function, and expression of HSP70 genes in B. rapa.

Identification and Characterization of the HSP70 Gene Family
In the first step of genome wide analysis of the HSP70 gene family to identify and characterization, the amino acid sequences of A. thaliana HSP70 family genes (AT3G12580.1) was obtained from the TAIR Arabidopsis genome database (TAIR-Home Page arabidopsis. org (accessed on 20 April 2022). This sequence was used as a query to retrieve the HSP70 genes in the Brassica database (http://brassicadb.cn/#/BLASTP/, accessed on 20 April 2022). Secondly, the HMMER 3.1 (http://www.hmmer.org, accessed on 20 April 2022) and BLASTP with the threshold e-value set to 1e −5 were performed using the Hidden Markov Model (HMM V 3.0) profiles of the HSP70 gene family (PF00012) were downloaded from the Pfam protein database (http://pfam.xfam.org/, accessed on 20 April 2022) and used as the inquiry, The default limitation of HMMER 3.1 was set to 0.01. Finally, total of 28 HSP70 were retrieved in the Brassica rapa genome. B.rapaHSP70 proteins were further characterized by determining the molecular weight, number of amino acids, isoelectric point, and grand average of hydropathicity (GRAVY) through the ProtParam tool (ExPASy-ProtParam tool, accessed on 20 April 2022). The protein sequences were uploaded to the online database Wolf PSORT for the prediction of subcellular localization (https:// wolfpsort.hgc.jp/, accessed on 20 April 2022). Furthermore, the protein conserved domain was identified using the NCBI conserved domain online server (www.ncbi.nlm.nih.gov, accessed on 20 April 2022) and motif analysis was performed by using MEME Suite (memesuite.org, accessed on 20 April 2022). The gene structure was predicted by the Gene Structure Display Server (GSDS) (http://gsds.cbi.pku.edu.cn/, accessed on 20 April 2022), using the CDS and genomic sequences of 28 selected B.rapaHSP70. Furthermore, prediction of the subcellular location pattern of each B.rapaHSP70 was carried out using the WoLF PSORT server (https://wolfpsort.hgc.jp/, accessed on 20 April 2022; [16]).

Phylogenetic Tree and Synteny Analysis of B.rapaHSP70 Family Proteins
The protein sequences of HSP70s from Arabidopsis, B. oleracea, and B. napus were retrieved to construct the phylogenetic tree using MEGA X (V 6.06) software (https:// www.megasoftware.net/, https://wolfpsort.hgc.jp/, accessed on 20 April 2022). The sequences were multiple aligned and employed using the neighbor-joining (NJ) method with 1000 bootstrap replicates [17]. Synteny relationships of B. rapa with B. oleracea, B. napus, and Arabidopsis were developed using MCScanX to obtain collinearity files that were used to build dual synteny plots in TBtools (https://github.com/CJ-Chen/TBtools, accessed on 20 April 2022).

B.rapaHSP70 Duplicated Gene Identification and Purity Selection Analysis
Duplicated genes were retrieved using the Blast, MCScanX, and Advance Circos features of TBtools [18]. For the purity selection pressure of duplicated genes, the Ks/Ka value was calculated using the Ka/Ks Calculator in TbTools, and TMY was calculated using previous studies in B. rapa [19,20]. Synteny analysis was performed for the correlation prediction among B. rapa and its closely related species through TBtools [21].

Cis-Element Analysis of B.rapaHSP70 Gene Promoters and Protein Structure
Putative cis-elements of the 28 B.rapaHSP70 genes were classified, and the 2000 bp upstream of the start codon were downloaded from the BRAD database (http://brassicadb. cn, accessed on 20 April 2022). The PlantCARE web-based tool (http://bioinformatics. psb.ugent.be/webtools/plantcare/html/, accessed on 20 April 2022) was used for further classification, and the results were presented using TBtools. The graphical construction of the 3-dimensional (3D) structure of B.rapaHSP70 proteins was done using the online web-based software PHYRE2 (PHYRE2 Protein Fold Recognition Server (ic.ac.uk), selecting the fold recognition method [22].

Targeted miRNA Prediction and Functional Analysis
CDS sequences of B.rapaHSP70 were used to determine the interaction of genes with microRNAs using the psRNATarget database (http://plantgrn.noble.org/psRNATarget, accessed on 20 April 2022), and the interactions were further prophesied by Cytoscape Software. Gene ontology (GO) annotation was performed to explore functional characteristics. The B.rapaHSP70 gene sequences were uploaded to the online eggNOG database (http://eggnog-mapper.embl.de/, accessed on 20 April 2022). Then, the GO annotation data were handled by using (http://brassicadb.cn, accessed on 20 April 2022) and graphical demonstration was given by using OmicShare (https://www.omicshare.com/, accessed on 20 April 2022).

Digital Expression Analysis of B.rapaHSP70
RNA-seq data were used to determine the expression of B.rapaHSP70 genes in different tissues, and a heat map was constructed using TbTools [23]. The data used for expression profiling were retrieved from the NCBI Gene Expression Omnibus (http://www.ncbi.nlm. nih.gov/geo/, accessed on 20 April 2022) under accession no. GSE43245. In addition, expression profiles of B.rapaHSP70 under cold and heat stress conditions were obtained by performing qRT-PCR.

Plant Material and Stress Conditions
Seeds of Chinese cabbage were grown in pots containing soil: vermiculite mixture (3:1) in a controlled chamber (24/18 • C, 60-65% relative humidity, 16/8 day/night) until they reached the five-leaf stage. Cold and heat stress conditions were applied to five-leaf seedlings. Seedlings were transferred to 4 • C for cold stress and 45 • C for heat stress under the same day/night and humidity conditions. Three biological replications were used in both stress conditions. Leaf samples were harvested after 0, 1, 3, 6, 12, and 24 h, placed in liquid nitrogen, and stored at −80 • C [24].

RNA Extraction and qRT-PCR Analysis
Leaf samples were used to extract the RNA, and qRT-PCR was performed with three replicates to check B.rapaHSP70 gene expression; actin was used as the internal reference gene in Chinese cabbage. All primers used in this experiment are listed in Table S4 The 2 −∆∆Ct method was used to check the relative expression levels, and the results were graphically represented [25].

Identification and Characteristics of HSP70 Gene Family Members in B. rapa
In this study, 28 genes were selected for further analysis using the A. thaliana query sequence (Gene I.D AT3G125801), and all selected genes were named B.rapaHSP70-1-B.rapaHSP70-28. These proteins were further analyzed for the HSP70 domain (PF00012) (Figure 1). Conserved motif along with phylogenetic tree demonstrate that all the B.rapaHSP70s are clustered into four groups (1, 2, 3 and 4), equal number and similar motif are present within each group except in group 4 ( Figure 1A). Detailed information about the basic physical characteristics of the proteins is given in Table 1. Within the 28 HSP70 proteins, 12 members were localized in the cytoplasmic membrane, 5 in the endoplasmic reticulum, 4 in the chloroplast, 4 in the mitochondria, and 3 in the nucleus (Table 1). Gene structure predicts that the genes localized in the nucleus consists of the higher number of intron and exon while the cytoplasmic genes consist of the lower number of intron and exon ( Figure 1B). The amino acid length of B.rapaHSP70 proteins was 636-894 aa. Among these, the shortest protein was B.rapaHSP70-15(636 aa), and the longest was B.rapaHSP70-28 (894 aa). The genes that translate into HSP70 proteins were distributed on 10 chromosomes, A01-A10. Among these, the maximum genes were distributed on chromosome A03 and the minimum on chromosomes A02, A04 and A05 ( Figure 2).

Evolutionary Relationship of B.rapaHSP70 Gene
The evolutionary relationships between B.rapaHSP70, B.nHSP70, B.olHSP70, and AtHSP70 were analyzed. Based on domains and a phylogenetic tree, 77 HSP70s were clustered into six major groups (A, B, C, D, and E) ( Figure 3). The results from the evolutionary relationship showed that group A consists of 11 B. rapa, among which 5 members belonged to Arabidopsis. Group B was the largest group, containing 10 members, among which 4 members belonged to B. rapa, 1 member to B. napus, 3 to B. oleracea, and 2 to Arabidopsis. A total of 17 HSP70 members were clustered into group C, 4 from B. rapa

B.rapaHSP70 Gene Duplication and Ka/Ks Analysis
Gene duplication is one of the most important characteristics of a plant genomic structure, and it can occur by independent mechanisms resulting in segmental or tandem duplications. Due to the importance of gene duplications in the evolution of gene families in plants, we analyzed the gene duplication of putative B.rapaHSP70 genes in the B. rapa genome. We detected duplicated gene couples among the 28 B.rapaHSP70 genes identified in B. rapa ( Figure  4A). Duplicated gene BraA08g021750.3C/BraA01g001130.3C, BraA03g047170.3C/BraA01g019900.3C, and BraA06g001280.3C /BraA04g004290.3C.The value of Ka/Ks can be used as an indicator of the selection pressure of a gene during evolution. All of the values were <1, suggesting that the B.rapaHSP70 genes were primarily evolved under the influence of purifying selection [26]. The divergence time million year ago (TMY) value was calculated based on T = Ks/2x, and the value of x was used according to [20] (Table 2).

B.rapaHSP70 Gene Duplication and Ka/Ks Analysis
Gene duplication is one of the most important characteristics of a plant genomic structure, and it can occur by independent mechanisms resulting in segmental or tandem duplications. Due to the importance of gene duplications in the evolution of gene families in plants, we analyzed the gene duplication of putative B.rapaHSP70 genes in the B. rapa genome. We detected duplicated gene couples among the 28 B.rapaHSP70 genes identified in B. rapa ( Figure 4A) The value of Ka/Ks can be used as an indicator of the selection pressure of a gene during evolution. All of the values were <1, suggesting that the B.rapaHSP70 genes were primarily evolved under the influence of purifying selection [26]. The divergence time million year ago (TMY) value was calculated based on T = Ks/2x, and the value of x was used according to [20] (Table 2).     Figure 4A), which shows that segmental duplication events occurred during evolution. Furthermore, a number of orthologs were identified among A. thaliana, B. rapa (AA), B. napus (AC), and B. oleracea (CC) using collinearity analysis. In addition, a comparative synteny analysis of HSP70 gene pairs among B. rapa, B. napus, B. oleracea, and A. thaliana was conducted ( Figure 4B). We selected these species on the base of highest morphotype similarity of B. rapa (such as Chinese cabbage) and B. oleracea (Cabbage) and B. napus (Mustard) was included due the presence of both sub genome (A and C). B.rapaHSP70 displayed the most collinearity with B. napus, followed by A. thaliana and B. oleracea. Results demonstrated that homologues genes from A genome of B. rapa are present in A as well as C genome of B.napus. Similar behavior was also observed in C genome of B. oleracea. These results suggest that, in addition to the whole genome duplication event, an independent duplication event also occurred during the evolution of these species.

Cis-Element Analysis in Promoter Regions of B.rapaHSP70 and Their Distribution
The sequence length of 2000 bp of all identified B.rapaHSP70 genes was retrieved and uploaded to the PlantCare database to examine the cis-elements in the promoter region. The graphical representation ( Figure 5A) illustrates a summary of the cis-acting elements. In total, 5 phytohormonal (ABA, abscisic acid; MeJA, methyl jasmonate, GA, gibberellin, and SA, salicylic acid) associated elements and 4 abiotic stress (temperature, drought, and light, and salinity) responsive elements were identified. Maximum elements were found for light response and auxin hormone response, and minimum elements were found for salicylic acid response in all promoters of the 28 B.rapaHSP70 ( Figure 5B and Table S1). Cells 2022, 11, x FOR PEER REVIEW 11 (A) (B)

Prediction of the 3D Structures of B.rapaHSP70s
With the advent of graphical visualization of genomic data, 3D protein structures also explain the properties, as well as facilitate the comparative studies, of proteins. All 28 B.rapaHSP70s structures elucidate that all the HSP70 members in B. rapa possess the uniform monomer structure with alpha (α) helices and beta (β) sheets, because proteins with similar structures often have similar functions ( Figure 6). Protein prediction also showed the existence of HSP70 domain (PF00012) in all the B.rapaHSP70s. The protein structures were modelled via c2khoA template at a confidence level of 100%. The overall fold of the B.rapaHSP70 domain remains similar to other known HSP70 with the typical (αβα) topology, thus preserving conserved signature fold of HSP70 family with central mixed β-sheet sandwiched between α-helices [27]. On the bases of protein structure guide design, binding ligand can be determined, which are essential to bind with the external or internal residues. In addition, the tertiary structure of proteins predicted physical connections, distinct functional roles, independent divergence sequences, and conserved biological properties in the B.rapaHSP70 protein family.

Prediction of the 3D Structures of B.rapaHSP70s
With the advent of graphical visualization of genomic data, 3D protein structures also explain the properties, as well as facilitate the comparative studies, of proteins. All 28 B.rapaHSP70s structures elucidate that all the HSP70 members in B.rapa possess the uniform monomer structure with alpha (α) helices and beta (β) sheets, because proteins with similar structures often have similar functions ( Figure 6). Protein prediction also showed the existence of HSP70 domain (PF00012) in all the B.rapaHSP70s. The protein structures were modelled via c2khoA template at a confidence level of 100%. The overall fold of the B.rapaHSP70 domain remains similar to other known HSP70 with the typical (αβα) topology, thus preserving conserved signature fold of HSP70 family with central mixed β-sheet sandwiched between α-helices [27]. On the bases of protein structure guide design, binding ligand can be determined, which are essential to bind with the external or internal residues. In addition, the tertiary structure of proteins predicted physical connections, distinct functional roles, independent divergence sequences, and conserved biological properties in the B.rapaHSP70 protein family. Figure 6. Predicted 3D models of B.rapaHSP70 proteins. Models have been generated by using the Phyre 2 server in intensive mode. Models were visualized by rainbow colour from N to C terminus.
Cells 2022, 11, x FOR PEER REVIEW 13 of 24 Figure 6. Predicted 3D models of B.rapaHSP70 proteins. Models have been generated by using the Phyre 2 server in intensive mode. Models were visualized by rainbow colour from N to C terminus.

Functional Annotation Analysis of B.rapaHSP70 Genes
Functional analysis using GO annotation and enrichment terms, such as molecular function (MF), cellular component (CC), and biological process (BP), were investigated to further explore the B.rapaHSP70 genes' functions. The annotation results for BP, MF, and CC exhibited several significantly enriched terms (Table S3). Some of the most enriched terms are highlighted in the subsequent section with a graphical illustration (Figure 8).  (Table S3).

Expression Profiling of B.rapaHSP70 Genes in Various Tissues
To illustrate the transcript levels of the B.rapaHSP70 genes, we examined 8 tissues and organs of B. rapa at various growth phases based on RNA-seq data of B. rapa ('Chiifu-401-42' variety; NCBI Gene Expression Omnibus; http://www.ncbi.nlm.nih.gov/geo/, accessed on 20 April 2022) under accession no. GSE43245. For instance, the expression patterns of most B.rapaHSP70 genes in the silique, callus, and flower were higher than those of other tissues (Figure 9). Three B.rapaHSP70 genes (B.rapaHSP70-10, B.rapaHSP70-11,  B.rapaHSP70-12) showed no transcript changes in any tissue/organ.  (Table S3).

B.rapaHSP70 Gene Expression in Response to Heat and Cold Stress
To examine the response of B.rapaHSP70s under cold and heat stress, we investigated the gene expression of B.rapaHSP70s in Chinese cabbage plants grown at 4 °C under cold stress and 45 °C under heat stress conditions. qRT-PCR was used to investigate the transcript levels of each B.rapaHsp70 gene with three biological repetitions and two technical  rapaHSP70-10, B.rapaHSP70-11, B.rapaHSP70-12, B.rapaHSP70-16, B.rapaHSP70-20, and B.rapaHSP70-26 showing no expression change. These findings suggest that these candidate genes may play diverse roles in regulating B. rapa growth processes (Figure 9).

B.rapaHSP70 Gene Expression in Response to Heat and Cold Stress
To examine the response of B.rapaHSP70s under cold and heat stress, we investigated the gene expression of B.rapaHSP70s in Chinese cabbage plants grown at 4 • C under cold stress and 45 • C under heat stress conditions. qRT-PCR was used to investigate the transcript levels of each B.rapaHsp70 gene with three biological repetitions and two technical repetitions. The primers used in this experiment are provided in Table S4. Generally, the relative expression level of the B.rapaHSP70 genes under all stress conditions fluctuated during the 24-h treatments.
Most of the B.rapaHSP70 genes were sensitive to heat stress, and most were highly regulated under high-temperature stress compared to cold stress. A few B.rapaHSP70s showed no significant expression differences after being treated under different time intervals. The expression levels of 16 B.rapaHSP70s (B.rapaHSP70-1, -2, -3, -4, -8, - 14, -15, -16, -19, -20, -22, -23, -24, -25, -27, and -28) were upregulated under high-temperature stress at different time intervals, and most were highly expressed after 24 h of heat stress. The expression levels of 16 B.rapaHSP70s (B.rapaHSP70-1, -2, -3, -4, -7, -8, - 13, -14, -17, -20, -21, -22, -23, -24, -26, and -28) were high at different time intervals under cold stress ( Figure 10 and Figure S5). qRT-PCR analysis showed that most of the B.rapaHSP70s expression were higher under heat stress as compared to cold stress. Interestingly, it is noticeable that B.rapaHSP70-3, B.rapaHSP70-14, and B.rapaHSP70-20 showed higher expression under both stress conditions. With the comparison of digital expression pattern results and our experimental analysis, there is significant difference in the expression of B.rapaHSP70 genes. The B.rapaHSP70 genes that are showing no expression changes in the normal condition such as B.rapaHSP70-10, B.rapaHSP70-11, B.rapaHSP70-12, and B.rapaHSP70-16, there expression is higher under the heat and cold shock conditions. All 28 B.rapaHSP70s exhibited a change in expression under stress conditions, these findings may help the breeders to improve the heat and cold tolerance in the Chinese cabbage. repetitions. The primers used in this experiment are provided in Table S4. Generally, the relative expression level of the B.rapaHSP70 genes under all stress conditions fluctuated during the 24-h treatments.
Most of the B.rapaHSP70 genes were sensitive to heat stress, and most were highly regulated under high-temperature stress compared to cold stress. A few B.rapaHSP70s showed no significant expression differences after being treated under different time intervals. The expression levels of 16 B.rapaHSP70s (B.rapaHSP70-1, -2, -3, -4, -8, -14, -15, - 16, -19, -20, -22, -23, -24, -25, -27, and -28) were upregulated under high-temperature stress at different time intervals, and most were highly expressed after 24 h of heat stress. The expression levels of 16 B.rapaHSP70s (B.rapaHSP70-1, -2, - 3, -4, -7, -8, -13, -14, -17, -20, -21, -22, -23, -24, -26, and -28) were high at different time intervals under cold stress (Figures 10  and S5). qRT-PCR analysis showed that most of the B.rapaHSP70s expression were higher under heat stress as compared to cold stress. Interestingly, it is noticeable that B.rapa-HSP70-3, B.rapaHSP70-14, and B.rapaHSP70-20 showed higher expression under both stress conditions. With the comparison of digital expression pattern results and our experimental analysis, there is significant difference in the expression of B.rapaHSP70 genes. The B.rapaHSP70 genes that are showing no expression changes in the normal condition such as B.rapaHSP70-10, B.rapaHSP70-11, B.rapaHSP70-12, and B.rapaHSP70-16, there expression is higher under the heat and cold shock conditions. All 28 B.rapaHSP70s exhibited a change in expression under stress conditions, these findings may help the breeders to improve the heat and cold tolerance in the Chinese cabbage.

Discussion
In this cutting-edge era of modern science, with the availability of open-source databases of whole genome sequences, HSP70 gene families have been explored in model plants and a diverse range of horticultural plant species, such as soybean. With the availabilities of the whole genome sequence of many plants, several HSP70 families have been identified, such as in A. thaliana and N. tabacum, as well as in different horticultural plants, such as B. napus, B. oleracea, and pumpkin [28][29][30][31][32]. However, this technique (genome-wide analysis) has not yet been employed in the HSP70 family in Brassica rapa.

Characterization of the B.rapaHSP70 Gene Family in B. rapa
The model plant A. thaliana belongs to the same family (Brassicaceae) as B. rapa. We identified 28 gene family members of HSP70 in B. rapa, which is higher than that found in Arabidopsis (11 AtHSP70) [8]. The identified genes were analyzed to explore their physical properties. There was no difference in the location of exon/intron boundaries or coding sequences. However, there were significant variances in the sequences and sizes of the introns in each of the 28 HSP70 genes.
Notably, genes within the same phylogenetic subgroup had similar motif compositions and exon/intron structures ( Figure 1 and Table 1). Genes without introns or those with few introns have been found to have enhanced plant expression levels in some studies [33]. In addition, the gene classifications of B.rapaHSP70 were confirmed by the correlation between intron numbers and motif arrangements in combination with phylogeny. A phylogenetic tree was created to show the evolutionary relationships between B. rapa, B. napus, B. oleracea, and Arabidopsis. Recently, research has examined the evolutionary associations among these species [32,34].
During the genome evolution process, gene duplications and chromosomal segments are major forces in the evolution of plant genome structure and content [35]. Tandem or segmental duplication events in which two or more neighboring genes duplicate on the same chromosome or on different chromosomes, respectively, can occur [36]. Based on the importance of duplication events, gene duplication events were explored to further understand the expansion mechanism of the B.rapaHSP70 family genes and we find all the duplicated B.rapaHSP70s are segmental duplicated.

Discussion
In this cutting-edge era of modern science, with the availability of open-source databases of whole genome sequences, HSP70 gene families have been explored in model plants and a diverse range of horticultural plant species, such as soybean. With the availabilities of the whole genome sequence of many plants, several HSP70 families have been identified, such as in A. thaliana and N. tabacum, as well as in different horticultural plants, such as B. napus, B. oleracea, and pumpkin [28][29][30][31][32]. However, this technique (genome-wide analysis) has not yet been employed in the HSP70 family in Brassica rapa.

Characterization of the B.rapaHSP70 Gene Family in B. rapa
The model plant A. thaliana belongs to the same family (Brassicaceae) as B. rapa. We identified 28 gene family members of HSP70 in B. rapa, which is higher than that found in Arabidopsis (11 AtHSP70) [8]. The identified genes were analyzed to explore their physical properties. There was no difference in the location of exon/intron boundaries or coding sequences. However, there were significant variances in the sequences and sizes of the introns in each of the 28 HSP70 genes.
Notably, genes within the same phylogenetic subgroup had similar motif compositions and exon/intron structures ( Figure 1 and Table 1). Genes without introns or those with few introns have been found to have enhanced plant expression levels in some studies [33]. In addition, the gene classifications of B.rapaHSP70 were confirmed by the correlation between intron numbers and motif arrangements in combination with phylogeny. A phylogenetic tree was created to show the evolutionary relationships between B. rapa, B. napus, B. oleracea, and Arabidopsis. Recently, research has examined the evolutionary associations among these species [32,34].
During the genome evolution process, gene duplications and chromosomal segments are major forces in the evolution of plant genome structure and content [35]. Tandem or segmental duplication events in which two or more neighboring genes duplicate on the same chromosome or on different chromosomes, respectively, can occur [36]. Based on the importance of duplication events, gene duplication events were explored to further understand the expansion mechanism of the B.rapaHSP70 family genes and we find all the duplicated B.rapaHSP70s are segmental duplicated.

Promoter Analysis and Gene Ontology
Previously, a total of 11 HSP70 genes from A. thaliana were identified and subjected to promoter analysis to investigate the presence of temperature-associated (heat shock element (HSE)) and dehydration-associated cis-elements (dehydration responsive element; DRE) [8]. Based on our findings, most HSP70 genes were predicted to be light responsive, presenting their effects when light or heat intensities were altered.
These stress-related cis elements revealed the involvement of HSP70 genes in B. rapa stress tolerance. Furthermore, GO enrichment analysis supports the findings of the HSP70 gene association with stress tolerance (Table S2). The molecular mechanism and functional investigation of stress-related HSP70 genes in Brassica crops should be a focus in the future to understand the underlying mechanisms.

miRNA: Key Players in the Regulation of Stress Responses
Several studies have revealed the importance and role of microRNAs (miRNAs) in post-transcriptional gene regulation [37]. In the past few years, numerous miRNAs have been identified in B. rapa that respond to different environmental stress conditions [13,38]. For instance, miR166 has been upregulate due to sunlight in maize crops and under heat stress in Chinese cabbage [39]. Similarly, miR393 has been reported to be associated with cold stress in tea plants (Camellia sinensis) [40]. These studies suggest that these Brassicaassociated miRNAs might play decisive roles against numerous stressors by modifying the transcript levels of HSP70 genes. In our study, we found that Bra-miR162, Bra-miR165, Bra-miR172, Bra-miR316, Bra-miR375, Bra-miR396, Bra-miR397, Bra-miR408, Bra-miR858, Bra-miR5714, and Bra-miR5717 families targeted the B.rapaHSP70 genes (Tables 3 and S3).

Differential Expression Analysis of B.rapaHSP70 Genes in Chinese Cabbage Leaves under Heat and Cold Stress
The expression patterns of genes are associated with their function [41]. Several studies have reported the use of qRT-PCR for investigating transcript levels, such as expression profiling of StHSP20 under heat stress. In another study, the results verified the involvement of B.rapaHSP70 proteins in thermotolerance [42]. Until now, few HSP70 genes have been identified and functionally described in B. rapa and other plant species. In this study, we laid a foundation for recognizing signaling controlled by HSP70 proteins under biotic and abiotic stress conditions.

Conclusions
In the current study, we performed a genome-wide analysis of the HSP70 gene family in B. rapa, and 28 putative members were identified. In silico analyses were performed, including gene structure, distribution, phylogenetic relationship, and syntenic studies, which helped to explore the evolutionary properties of the HSP70 gene family in B. rapa. In addition, targeted miRNAs, promotor cis-acting regulatory elements, and gene ontology (GO) were executed. The finding of all these analysis elucidating that B.rapaHSP70s are targeted by 34 different families of microRNAs, these genes are highly responsive to light, temperature and phytohormones, and B.rapaHSP70s are involved in the Biological, cellular and molecular functioning of B. rapa, respectively. Three dimensional protein structural knowledge can help improve crop properties, such as improving stress resistance and biomass yield. Furthermore, quantitative real-time PCR results indicated that HSP70 genes are strongly involved in heat and cold stress responses in B. rapa plants. To a large extent, all analyses performed on the HSP70 gene family will lay the foundation for further studies of molecular and physiological functions in Brassica crops.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cells11152316/s1, Table S1; Cis-Element Analysis of B.rapaHSP70 Gene Promoters. Putative cis-elements of the 28 B.rapaHSP70 genes identified via PlantCARE webbased tool, Table S2; Genome-Wide Analysis of miRNA-Associated B.rapaHSP70 Genes. Detailed information of the miRNA-targeted sites determined by psRNATarget database, Table S3; Functional Annotation Analysis of B.rapaHSP70 Genes. Detailed information of the GO annotation data retrived from Brassica database, Table S4; Primers used in qRT-PCR. Expression patteren analysis in Response to Heat and Cold Stress, Table S5

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request. However, most of the data is shown in Supplementary files.