Genome-wide characterization and expression profiling of PDI family gene reveals function as abiotic and biotic stress tolerance in Chinese cabbage (Brassica rapa ssp. pekinensis)

Background Protein disulfide isomerase (PDI) and PDI-like proteins contain thioredoxin domains that catalyze protein disulfide bond, inhibit aggregation of misfolded proteins, and function in isomerization during protein folding in endoplasmic reticulum and responses during abiotic stresses.Chinese cabbage is widely recognized as an economically important, nutritious vegetable, but its yield is severely hampered by various biotic and abiotic stresses. Because of, it is prime need to identify those genes whose are responsible for biotic and abiotic stress tolerance. PDI family genes are among of them. Results We have identified 32 PDI genes from the Br135K microarray dataset, NCBI and BRAD database, and in silico characterized their sequences. Expression profiling of those genes was performed using cDNA of plant samples imposed to abiotic stresses; cold, salt, drought and ABA (Abscisic Acid) and biotic stress; Fusarium oxysporum f. sp. conglutinans infection. The Chinese cabbage PDI genes were clustered in eleven groups in phylogeny. Among them, 15 PDI genes were ubiquitously expressed in various organs, while 24 PDI genes were up-regulated under salt and drought stress. By contrast, cold and ABA stress responsive gene number were ten and nine, respectively. In case of F. oxysporum f. sp. conglutinans infection 14 BrPDI genes were highly up-regulated. Interestingly, BrPDI1–1 gene was identified as putative candidate against abiotic (salt and drought) and biotic stresses, BrPDI5–2 gene for ABA stress, and BrPDI1–4, 6–1 and 9–2 were putative candidate genes for both cold and chilling injury stresses. Conclusions Our findings help to elucidate the involvement of PDI genes in stress responses, and they lay the foundation for functional genomics in future studies and molecular breeding of Brassica rapa crops. The stress-responsive PDI genes could be potential resources for molecular breeding of Brassica crops resistant to biotic and abiotic stresses. Electronic supplementary material The online version of this article (10.1186/s12864-017-4277-2) contains supplementary material, which is available to authorized users.


Background
Protein disulfide isomerases (PDIs) are enzymes found primarily in the endoplasmic reticulum (ER) in eukaryotes play a vital role in protein folding. Proteins, immediately after biosynthesis must be folded into correct threedimensional shape for proper functioning. Misfolded proteins are often nonfunctional by producing aggregates and interferes cellular functions. PDIs can bind into misfolded or unfolded proteins preventing to produce such aggregates [1]. Venetianer and Straub, [2] firstly identified PDIs as protein-folding catalysts and demonstrated their catastrophic consequences of a defective protein folding process. Plants utilize multiple mechanisms to fold proteins correctly for proper functioning. PDIs are one of the mechanisms, which work through formation and breakage of connections between cysteine residues by producing disulfide bonds, which was described by Anfinsen [3]. These bonds stabilize the folded protein and provide correct structure to perform its particular function.
PDIs contain four thioredoxin (TRX)-like domains, two of which contain a catalytic site for disulfide bond formation. The reduced (dithiol) forms of PDIs catalyze the reduction of impaired thiol residues of a particular substrate, acting as an isomerase [4]. In human, PDIs usually contain four TRX-like domains (a, b, b' and a'), a linker (x) and a C-terminal extension domain (c) [5]. The 'a' and 'a' domains are TRX domains, containing an active Cys-Gly-His-Cys motif joined to α-helices and βstrands (β-α-β-α-β-α-β-β-α) that is essential for polypeptide redox and isomerization [6]. Although the 'b' and 'b' domains share some structural identity with TRX domain, they do not possess specific active motif site [6]. The α-helices are particularly important for DNA binding motifs, including helix-turn-helix motifs, leucine zipper motifs and zinc finger motifs. The helix-turn-helix (HTH) motif is a major structural motif capable of DNA binding. This motif is composed of two α-helices joined by a short strand of amino acids found in many regulatory proteins for gene expression [7]. The C-terminal region contains a C-domain rich in acidic residues typical as calcium binding proteins and ends with an ER retention signal composed of four amino acids, generally KDEL/GKNF/VASS [8].
PDI genes have been identified in many higher plants, such as 21 PDI genes in Arabidopsis thaliana, 12 in Oryza sativa, 12 in Zea mays, 10 in Vitis vinifera, 9 in Triticum aestivum and so on [9]. In eukaryotes, PDI family proteins are divided into eleven groups based on phylogenetic analysis. Proteins in groups' I-V contain two thioredoxin-like active domains, whereas those in groups VI-XI contain a single thioredoxin-like active domain. PDI proteins in higher plants are involved in signal transduction pathways and in transcriptional complexes that regulate the responses of genes to environmental stimuli. PDI proteins ERP57, PDIp, P5, ERP72, PDIR, and PDI-D act as redox catalysts and isomerases and exhibit differential functions, such as peptide binding, cell adhesion, and chaperone activities [10]. These proteins play key role as storage protein and plasma membrane maturation [9] GmPDIL-3a and GmPDIL-3b are highly expressed during seed maturation suggesting that they are involved in folding or accumulation of storage proteins in soybean [11]. PDIL2-1 of A. thaliana is directly involved in ovule structure and embryo sac development and determining proper direction of pollen tube growth [12]. The PDI-like protein RB60 plays important roles providing redox potential to regulate photosynthesis [13]. Most wheat PDI and PDIL genes are expressed during endosperm development, indicating their association with storage protein biosynthesis and deposition, which is directly related to gluten quality [14]. BdPDIL1-1 is significantly up-regulated under abiotic (drought, salt, ABA and H 2 O 2 ) stress, suggesting their involvement in multiple stress responses [15]. All maize PDI genes are highly responsive to cold, salt, dehydration, and ABA stress [16]. During drought, heat and cold stress, TaPDI1, TaPDI2 and TaPDI3 are highly up-regulated in roots and other tissues [17]. There are many family genes involved in tolerance to abiotic and biotic stresses in plants; however we have selected PDI because of its abundance of ubiquitous sulfydryl oxidoreductase protein, which is an important cellular protein with multiple biological functions of all eukaryotic cells. This sulfydryl oxidoreductase protein of PDI displays versatile redox behavior, which also interact with other proteins and assumed its potential role against various diseases and abiotic stresses [18]. In addition, PDI is an important redox proteins regulating reactive oxygen species (ROS) production in the cells and alter redox status of cells to activate defense mechanism [19]. In this study, we make effort to identify the members of PDI gene family in Brassica rapa through genome wide exploration by using different bioinformatic tools. In addition, the putative candidates of PDI genes responsive to abiotic and biotic stresses are predicted through expression profiling using stress induced Chinese cabbage materials.

Results
Identification of PDI genes in B. rapa A total of 32 BrPDI genes were identified using SWIS-SPROT of the B. rapa genomic database (http://brassicadb.org/brad/searchAll.php) (BRAD; [20]); using key word "PDI", NCBI, and annotations of microarray Br135K dataset from cold-treated B. rapa (Chiifu & Kenshin), removing any duplicates. A BLAST search was performed using Arabidopsis PDI sequences as the query, and picked the encoded protein and genomic sequences of PDI genes from BRAD database for 32 BrPDI genes. A high degree of similarity of these 32 BrPDI protein sequences was also picked for other plant species. Isoelectric points, molecular weights and residual size (ranged from 144 to 596 aa) of putative 32 BrPDI proteins are presented in Table 1.

Phylogenetic analysis and domain location
A phylogenetic tree was constructed using 76 full-length protein sequences of PDI and PDI-like genes, which included 32 from Chinese cabbage, 11 from Brachypodium distachyon, 12 from maize, and 21 from Arabidopsis to investigate evolutionary relationship. Eleven phylogenetic groups were denoted among the considered PDI genes from different plant species ( Fig. 1) with distinct 4 clades. Clade 1 consisted phylogenetic groups I, II, III and VII, whose members contain two active thioredoxin domains, whereas members of group VII encoded proteins with a single N-terminal active domain (Fig. 2). Clade 2 contained phylogenetic group IV, V and VI, among them members of group IV and V possessing two active thioredoxin domains in tandem at their Ntermini. Group members of VI may have lost one of the two active thioredoxin domains due to shared structural features of a common progenitor. Clade 3 contained phylogenetic groups VIII, which included genes encoding proteins with a single active thioredoxin domain. The members of groups IX, X and XI formed clade 4; encoded proteins with an active thioredoxin domain. Phylogenetic and domain analyses revealed that these PDI family genes are divergent in plant. Most BrPDI proteins are contained an N-terminal signal peptide and a C-terminal KDEL signal responsible for translocation and ER retention, respectively. In the domain structures, 'a' and 'a' are homologous to thioredoxin and contain a -CXXC-active site for isomerase and redox activities. By contrast, the 'b' and 'b' domains have no homology to thioredoxin and lack of -CXXC-catalytic site. However, the secondary structures of four (a, a' , b, b') domains are similar to thioredoxin rather than the active catalytic site. The thioredoxin-like domain comprises βαβαβαββα motifs forming four layers of β-sheets that are sandwiched three layers of α-helixes (Fig. 3).
We performed multiple sequence alignment of 'a'and 'a'-type domain sequences of BrPDI proteins to analyze the features of them. These domain sequences contained five β-sheets, four α-helices and the -CGHC-motif of the active site. They also contained conserved arginine, glutamic acid, proline, and lysine residues (Fig. 3). All BrPDI members 'a'-type domain with conserved arginine residue, critical for their catalytic function, except BrPDI3-1 and BrPDI3-2 which contained phenylalanine or leucine residues instead of arginine residue. In addition, BrPDI9-1 contained tryptophan and BrPDI10-1, BrPDI10-2 and BrPDI10-3 contained lysine in the place of arginine residues. BrPDI2-3 accomplished with RGHC non-characteristic active sites; but BrPDI3-1 and BrPDI3-2 with CARS. Whereas, BrPDI8-1, BrPDI8-2, BrPDI8-5, and BrPDI8-6 comprised with CYWS; and BrPDI10-1, BrPDI10-2, BrPDI10-3, and BrPDI10-4 contained CPFS non-characteristic active sites (Fig. 2) which may affect their redox potential and lose their function. Chromosomal distribution of BrPDI genes Thirty two BrPDI genes are located on different chromosomes not less than one on each except chromosome 04. The highest number (six) of BrPDI genes was identified on chromosome A05 (19.35%), while chromosome A01, A02, A07, and A10 contained two BrPDI genes (6.45%) for each, and chromosome A04 has no BrPDI genes (0%). Only one gene (BrPDI11-3) is located in scaffold ( Fig. 4a and b). None of the BrPDI gene clusters was detected in Chinese cabbage. The sequence alignment of BrPDI proteins showed higher similarity within groups. BrPDI proteins showed ≥68% similarity within groups, except group10 and group 2. By contrast, similarity between groups was ≤65%. Seventeen pairs of BrPDI proteins showed 80% similarity indicating duplications are predominant for those genes ( Table 2). domain are shown. The thioredoxin-like catalytic domains with two active sites (shown in detail in the boxes) are also shown. Numbers above indicate domain boundaries (aa), and numbers on the right indicate ORF (aa). Domains a o (light gray) and a' (gray) are homologous to TRX and contain the catalytic CxxC motif (red). Domains b (yellow) and b' (light yellow) also exhibit a TRX fold, but they do not share high sequence similarity with each other or with domains a or a'. The C-terminal extension (red) contains a (K/H)DEL retention signal for the ER Three fractionated subgenomes, like least fractionated (LF), medium fractionated (MF1), and most fractionated (MF2) subgenomes are found in B. rapa genome. Notably, 32 BrPDI genes are distributed onto ten chromosomes with fractionation into three subgenomes, of which are 15 LF (46.87%), 8 MF1 (25%), and MF2 9 (28.13%; Fig. 4c and Table 1). In addition, 32 BrPDI genes are distributed in different block during evolution. Among them, 7 genes (22.58%) belong to AK8 block, followed by 19.45% of BrPDI genes into AK1, AK3 blocks; while only 6.45% of BrPDI genes were allocated to AK2, AK4, AK5 and AK6 (Fig. 4d). We have also found a total of 17 segmental duplicated PDI genes pairs in B. rapa genome ( Fig. 5 and Table 3). Furthermore, the substitution rate of non-synonymous (Ka) and synonymous (Ks) ratios (Table 3) were calculated to assess the selection pressures among duplicated BrPDI gene pairs. In these analysis, Ka/Ks ratios were identified as lower than 1, 1 and above 1 indicating negative or purifying selection, neutral selection and positive selection, respectively. Nine out of seventeen BrPDI duplicated gene pairs, had Ka/Ks value below 1 (purifying selection) and rest of them had Ka/Ks value above 1(positive selection). Moreover, the estimated divergence time of BrPDI genes showed that duplication event started at 29.74 million years ago (MYA) and continued up to 1.62 MYA (Table 3).

Motif composition and intron-exon analysis
The intron-exon structures were almost consistent among the groups of BrPDI genes. Group I contained nine to ten exons; whereas, members of group V contained nine exons (Additional file 1: Figure S1). All members of group II possessed 10 to 11 exons, while members of group IV and IX contained 11 exons, and the member of group VIII composed with the highest number (15) of exons. Member of groups X and XI possessed 3-4 exons. Finally, member of groups VI and VII contained four and five exons, respectively. Conserved motifs among BrPDI proteins were investigated using the MEME motif search tool. Motifs 1 and 2 contained-CXXC-catalytic sites, which are necessary for isomerase and redox activity. Motifs 3 was present in all groups, but motifs 2 was absent in groups VI, VII, VIII, IX and X. Motif 1 and 5 were absent in group IV and XI, respectively. Relatively less conserved motif 6 was found only in groups I, II, III Fig. 3 Multiple sequence alignment of the a-type domains of B. rapa PDI proteins and a typical human PDI. Residues highlighted in gray and light gray share 100% and >50% identity, respectively. Elements of the secondary structure are specified by blue colour bars (α-helices) and red colour bars (β-sheets). Red arrows indicate the two buried charged residues in the vicinity of the active site, orange arrow indicates the conserved arginine (R) and green arrow indicates the cis prolines (P) near each active site. Active-site residues within a domain are pink colour boxed and VII. Motifs 4, 7, 8 and 9 were unique to group VIII and motif 10 was present in group VIII and XI only (Additional file 1: Figure S2).

Microsyntenic analysis
A microsynteny map was constructed to find out orthologous gene pairs of PDI genes among B. rapa, B.distachyon and A. thaliana for exploration of evolutionary history and relationships among the genomes (Fig. 5). We found 22 orthologous gene pairs between B. rapa and A. thaliana; whereas, 13 orthologous gene pairs were identified between B. rapa and B. distachyon (Fig. 5). We found 17 duplicated PDI gene pairs in B. rapa genome. All of the segmental duplications are denoted with an orange line in fig. 5, no tandem duplications were found in PDI genes in B. rapa.

Microarray data analysis
We used our previously published microarray data set to analyze the expression patterns of 32 BrPDI genes using two contrasting inbred lines 'Chiifu' and 'Kenshin' , treated with cold and freezing stress (4°C, 0°C, −2°C and −4°C) for 2 h [20]. All of the BrPDI genes were differentially expressed in response to cold or freezing stress in two lines (Fig. 6). In Chiifu, most of the BrPDI genes were highly expressed upon exposure to cold and freezing temperatures. While, BrPDI2-1, 4-1, 8-2, 8-6, 9-1 and 10-1 genes showed higher expression during cold stresses in 'Kenshin' compared to 'Chiifu'.

Organ-specific gene expression
We have used gene-specific primers of 32 BrPDIs for RT-PCR to check the expression patterns of BrPDIs in different organs (root, stem, leaf, flower buds, sepals, petals, stamens, and pistils) in Chinese cabbage line SUN-3061. Semi-quantitative RT-PCR analysis revealed that among the 32 genes, 15 were ubiquitously expressed in all organs (Fig. 7). Three genes (BrPDI 1-2, 9-1 and 11-4) were slightly expressed in all of the tested organs. BrPDI 7-1 was highly expressed in floral parts but slightly expressed in roots, stems, leaves and flower buds. On the other hands, BrPDI 1-2, 1-3, 1-4, 2-3 and 8-1 were totally absent in roots, stems and leaves but only expressed in flower buds and floral parts. In addition, BrPDI 2-1 was weakly expressed in roots, leaves and pistil; while BrPDI3-1 was expressed slightly in roots, stems, and leaves. BrPDI8-5 was expressed in leaves, flower buds, and pistils; the BrPDI 8-6 was expressed in stamens only. BrPDI10-4 was appeared with very low signal in leaves, sepals and pistils. BrPDI 11-1 and BrPDI1 11-2 were highly expressed in vegetative organs but weekly expressed in reproductive organs (Fig. 7).

Expression profiling of the BrPDI genes during chilling injury treatment
Chilling injury experiment with cold tolerance 'Chiifu' and cold sensitive 'Kenshin' was conducted to get the deep insight of the expression of the target cold tolerance gene as an alternative of functional analysis of the predicted genes. We exposed the 'Chiifu' and 'Kenshin' seedlings in cold injury treatment at 0°C for 72 h until complete chilling injury of the susceptible (Kenshin). First remarkable chilling injury was recognized in 'Kenshin' seedlings at 24 h time point and gradually progress with the enhancement of time courses, 'Kenshin' seedlings were completely injured at 72 h time course. Whereas, no chilling injury was observed in 'Chiifu' seedlings up to 72 h time point, but plant growth become stunted (Fig. 9). In 'Chiifu' ,

Analysis of biotic stress-related gene expression
To elucidate the expression patterns of BrPDI genes in response to Fusarium oxysporum f. sp. conglutinans infection, we have collected leaf samples from infected and mock-treated Chinese cabbage line 'Chiifu' for performing qPCR. This fungus causes wilt and root rot diseases in Brassica crops. The BrPDI genes exhibited distinct expression patterns in response to Fusarium infection. Out of 32 genes, 14 showed high expression at 6 h time point (Fig. 11). The expression levels of BrPDI

Discussion
In genome-wide exploration study of Chinese cabbage PDI genes; we have identified 32 PDI genes identical to 21 PDI genes of Arabidopsis. The Chinese cabbage genome has gone under two duplication events since its divergence from Arabidopsis [21]. The evolutionary relationship between Arabidopsis and Chinese cabbage is also supported by our results. We analyzed 32 BrPDI genes expression pattern based on the microarray data set in two contrasting Chinese cabbage inbred lines, 'Chiifu' and 'Kenshin' exposed to cold and freezing stress (4°C, 0°C, −2°C, and −4°C) [20]. From an evolutionary point of view, gene duplications, tandem and segmental duplications can increase the gene numbers of a particular gene family [22]. Distribution of BrPDI genes in the B. rapa genome are also affected by segmental duplication, tandem duplication, polyploidization etc. during evolution [23,24]. In addition, genome triplication events of B. rapa may also be played important role in the expansion of PDI gene family. We found 17 pairs of segmental duplicated PDI genes in Chinese cabbage genome ( Fig. 5 and Table 3), this result indicates that the expansion of PDI gene family in B. rapa genome depicted on segmental gene duplication.
We also investigated evolutionary history of BrPDI gene family and calculated the Ka (the number of synonymous substitutions per synonymous site), Ks (the number of nonsynonymous substitutions per nonsynonymous site) and Ka/Ks ratios of segmental duplicated gene pairs. Nine gene pairs had Ka/Ks < 1 (purifying selection) and eight gene pairs had Ka/Ks >1(positive selection; Table 3) indicating purifying selection and positive selection together accelerated evolution and functional divergence of BrPDI genes in B. rapa genome. Based on the synonymous substitution rate, we calculated divergence time of the BrPDI genes and found that duplication started on 29.74 MYA and stop at 1.62 MYA (Table 3) which indicate PDI genes divergence take place after the genome triplication events in B. rapa. Genome triplication event of B. rapa has occurred from the ancestor of A. thaliana between five to nine million years ago (MYA; [25]). Besides this, the exons and introns distribution of BrPDI genes are also consistent with the number of exons and introns of PDI genes in Arabidopsis, which indicate their close evolutionary relationship. In addition with, by drawing microsyntenic map, we conclude that BrPDI genes are strongly related to those of B. distachyon and A. thaliana.
The Chinese cabbage PDI genes appeared as functionally differentiated for short time, and keeps as functional genes in the genome and maintaining complex functions in organ development. The features and functions of PDI family genes have been extensively studied in Triticum aestivum [14] and Brachypodium distachyon [15]. Our identified and characterized PDI genes in Chinese cabbage help to increase the understanding of role of this family gene. Most of the BrPDI genes were predominantly expressed in all organs suggested their functions in regulation of plant growth and development. Based on our organ-specific expression analysis, all BrPDI genes are expressed at least in one organ at variable level. BrPDI 2-1 and BrPDI 9-1 were highly expressed in root, which is consistent with the previous findings [14,15]. In addition to their involvement in plant growth and development, PDI genes are also responsive to adverse environmental conditions [15]. We have found BrPDI 1-2, BrPDI 1-3, BrPDI1-4, BrPDI 2-3 BrPDI 8-1 and BrPDI 8-6 were expressed in flower buds and floral parts, suggesting their possible role in flower development.
The expression patterns of 32 BrPDI genes were analyzed using a whole-genome microarray dataset of two inbred lines of B. rapa, 'Chiifu' and 'Kenshin' imposed to various temperature treatment [26]. Thereof, we have selected BrPDI family genes for analyzing their differential expression patterns compared to control by constructing heat-map (Fig. 6). Our results suggest that BrPDI genes play a vital role against abiotic stress responses in Chinese cabbage. In addition, our main objectives were to identify putative candidate PDI genes related to stress response during plant growth and development. It is established fact that environmental stresses like, cold, salt, and drought are severely affected crop production and make threaten in food security worldwide. Therefore, application of artificial stress treatment helps to know, how plants adapt to stresses via molecular, morphological, and physiological mechanisms. There is intense interest in identifying stress-responsive genes and elucidate them to develop stress-tolerant crop cultivars. Relative expression of PDI family genes were estimated from the sample of the plants subjected to abiotic stresses (cold, salt and drought), hormone (ABA), chilling injury, and biotic (F. oxysporum f.sp.conglutinans) stress treatment over different time periods (Fig. 8a-d,   9~11) using qPCR data. In plants, PDIs are localized to the ER and other cellular compartments, such as nucleus and chloroplasts [27,28]. Proteins in the ER of plants can lead to unfolding or misfolding due to complex physiological processes during different environmental stresses [29]. Plants possess multiple mechanisms for ensuring correct folding of proteins. PDI was the first reported catalyst of protein folding (2). In the current study, BrPDI1-3, BrPDI3-1, BrPDI6-1, BrPDI9-2, BrPDI10-2, and BrPDI11-1 were highly expressed in cold-tolerant 'Chiifu' than in cold-susceptible 'Kenshin' during cold stress treatment. PDI genes that are highly expressed under cold stress may encode proteins catalyze for breakage of disulfide bonds in misfolded proteins, restoring their proper functioning through correction of folding patterns. Lu & David [30] reported six Arabidopsis genes; those were up-regulated due to chemicals such as dithiothreitol (DTT), β-mercaptoethanol and tunicamycin induction. In current study, we have found BrPDI1-1, BrPDI4-1, BrPDI4-2, and BrPDI5-2 were highly expressed under cold stress in coldsusceptible 'Kenshin' compared to cold-tolerant 'Chiifu' , suggesting that these four genes might be related to the cold-susceptibility in 'Kenshin'. However, BrPDI2-1, BrPDI2-2, BrPDI5-2, BrPDI8-1, BrPDI8-6, BrPDI9-1, BrPDI10-4, and BrPDI11-4 were not expressed or even down-regulated in 'Chiifu' under cold treatment confirmed that these genes contain putative cold-inducible LTR cis-acting elements. Notably, these genes also contain a heat-inducible HSE cis-acting element. Hence, both cold-and heat-inducible cis-acting elements are presented in the promoter regions of the same genes suggesting that the latter cis-acting element may downregulate these genes under cold stress. Two contrasting cold tolerant Chinese cabbage lines were used in this study to predict the putative cold susceptible and cold tolerant BrPDI gene(s) of this crop. Four BrPDI genes (BrPDI1-4, BrPDI6-1, BrPDI9-2, and BrPDI11-1) exhibited higher expression in 'Chiifu' compared to 'Kenshin' at the level of chilling injury. Among these four genes BrPDI6-1, BrPDI9-2, and BrPDI11-1 genes were common in the higher expression in cold stress and chilling injury level treatments revealed their active involvement to overcome the cold and chilling injury at different time points of the cold tolerant line 'Chiifu'. Therefore, from the in depth expression data, it is clearly evident that BrPDI6-1, and BrPDI9-2 genes are the putative candidate in Chinese cabbage for overcoming cold and even chilling injury temperature (0°C).
With few exceptions, BrPDI genes are highly expressed in Chinese cabbage under salt and drought stress, among them BrPDI1-1 to BrPDI5-2 genes encode the proteins which containing two active sites and a Cterminal ER retention signal composed of four amino acids (KDEL/ VASS/ MADD), might be involved in enhancing protein-folding in the ER. BrPDI1-1, BrPDI1-3, BrPDI1-5, BrPDI3-1, BrPDI4-2, BrPDI5-1, BrPDI5-2, BrPDI8-5, and BrPDI8-6 were highly expressed under ABA stress compared to control. Mittler et al. [31] and Sakamoto et al. [32] observed rapid accumulation of reactive oxygen species (ROS) in plant cells subjected to environmental stresses such as salt, drought, and ABA. Indeed, during environmental stress (e.g., drought, salt, ABA, UV, and heat exposure) ROS levels can increase dramatically, which may significantly damage cell structure. PDI proteins play a significant role in thioredoxinbased redox pathway, which comprises part of the antioxidative defense system [33]. Abiotic stress usually leads to protein unfolding, misfolding, and aggregation; which represent a common threat to living cell [34]. Efficient protein repair systems and/or protein folding stability helps plant to survive in adverse environmental conditions [35]. Most BrPDI genes were differentially expressed in response to salt, drought, and ABA treatment, which is consistent with the findings of Zhu et al. [15]. Those genes are highly expressed in response to salt, drought, and ABA might function in the repair of misfolded proteins, thereby playing a role in plant resistance to the stresses. The up-regulated genes contain the respective stress-responsive cis-elements in their promoter regions (Additional file 2: Table S1).
PDI genes have been previously shown to be contributed against powdery mildew resistance in common wheat [36]. In this study, we have analyzed expression profiles of these genes in response to Fusarium oxysporum f.sp. conglutinans infection. Moreover, 14 genes were up-regulated compare to mock treated one due to infection of fungal pathogen at 6 h time point, whose genes containing a disease resistance Box-W1 cis-acting element (Fig. 11, Additional file 2: Table S1). By contrast, BrPDI 6-1~BrPDI 11-4 do not contain the Box-W1 element, those exhibited almost no expression due to fungal infection. Thus, the highly expressed BrPDI genes might play a vital role in response against fungal pathogen. Therefore, this might be the putative candidate genes for developing Fusarium oxysporum f.sp. conglutinans resistant Chinese cabbage cultivars through marker assisted back crossing (MAB) or gene technology. Besides, the up-regulated BrPDI genes might play differential roles in signal transduction pathways and/or cooperate with other genes to form networks for defending plants against adverse environmental conditions.

Conclusion
In conclusion, this is the first report of genome-wide characterization of PDIs in Chinese cabbage. We identified 32 BrPDI genes in the Chinese cabbage genome and characterized them based on motif distribution, protein structure, classification, number of introns -exons, sequence similarities, and expression patterns in response to abiotic (cold, salt, and drought) stresses, ABA treatment and biotic stress (F. oxysporum f.sp. conglutinans inoculation). We have also recognized BrPDI genes may play roles in biotic and abiotic stress responses. Several BrPDI genes (BrPDI1-1, BrPDI1-4, BrPDI1-5, BrPDI3-1, BrPDI3-2, BrPDI4-1, BrPDI4-2, BrPDI5-1, BrPDI5-2, BrPDI6-1, BrPDI7-1, BrPDI8-2, BrPDI9-1, BrPDI9-2, and BrPDI11-1) might function in responses to multiple stresses. The identified stress-induced BrPDI genes help to elucidate the complex regulatory network underlying stress resistance mechanisms. In addition, these genes represent a useful resource to molecular breeders for marker assisted back crossing (MAB) and/ or engineering transgenic plants in future research with increased resistance to biotic and abiotic stresses.

Identification and sequence analysis of BrPDI genes
Chinese cabbage PDI family members were identified using the SWISSPROT tool of the B. rapa database (http://brassicadb.org/brad/searchAll.php; 20) and NCBI using the keyword "PDI". We also searched these genes from a microarray Br135K dataset of two contrasting Chinese cabbage inbred lines, Chiifu and Kenshin exposed to low-temperature stress. The Arabidopsis PDI sequences were used as the query to perform a BLAST search setting a cutoff e-value of <10 −10 . The CDS (coding DNA sequences) and protein sequences of the identified PDIs were obtained from the B. rapa genomic database. The PDI protein sequences were further analyzed for the presence of a thioredoxin domain using the web tool "SMART" (http://smart.embl-heidelberg.de/ [37]) and NCBI (https://www.ncbi.nlm.nih.gov/cdd). Additionally, the primary structures of the genes (protein length, molecular weight and isoelectric point) were analyzed using ExPasy (http://web.expasy.org/compute_pi/). ORFs were identified using ORF finder at NCBI (https://www.ncbi.nlm.nih.gov/orffinder/). Sequence alignment of PDI proteins were carried out using CLUSTAL Omega (http://www.ebi.ac.uk/ Tools/msa/clustalo/). Sub-genome fractionation and positional information for all idetified BrPDI genes on 10 chromosomes of B. rapa were retrieved from the Brassica database, and draft maps of the locations of the BrPDI genes were constructed using Map Chart version 2.30 (http://www.wageningenur.nl/en/show/ Mapchart.htm). The Gene Structure Display Server (GSDS) web tool (http://gsds.cbi.pku.edu.cn/) was used to determine the number of introns and exons by aligning CDS and genomic sequences of the PDI genes [38]. Putative cis-acting regulatory elements in the BrPDI genes were predicted in regions of approximately1500bp upstream of the translation initiation site [ATG] using the PlantCARE (http://bioinformatics.psb.ugent.be/ webtools/plantcare/html/).

Phylogenetic and conserved motif analysis of Chinese cabbage PDI genes
The predicted BrPDI protein sequences were aligned with those of Arabidopsis, Brachypodium distachyon and Zea mays and Chinese cabbage using Clustal Omega, and a phylogenetic trees was constructed based on the condensed alignment in MEGA6.06 using the Neighbor-Joining (NJ) algorithm [39], with the parameters set at 1000 replications for bootstrap values and complete deletion mode to analyze tree topology and reliability. The genes which were used in phylogenetic analysis, name and accession number provided in Additional file 2: Table S2. Motif analysis of proteins was performed using MEME (Expectation Maximization for Motif Elicitation v4.10.1) [40] with the following parameters: (1) number of motifs = 10, (2) Motif width ≥ 6 and ≤50.

Chromosome localization and gene duplications
The physical locations of BrPDI genes were collected from B. rapa genomic database and the positions of the BrPDI genes were drafted to ten B. rapa chromosomes by map chart program. Duplicated BrPDI genes were identified using BLAST searched against themselves, notably identity and query coverage were >80% of those particular genes [41]. Tandem duplicated genes were marked as an array of two or more homologous genes within a distance of 100 kb. The synonymous rate (Ks), nonsynonymous rate (Ka), and evolutionary constraint (Ka/Ks) were calculated between the duplicated PDI gene pairs (Table 3), using the method of Nei & Gojobori [42] by Mega 6.0 software. The value of Ka/Ks ratio like >1, <1 and =1 are depicted for positive selection, purifying selection and neutral selection, respectively [43]. The divergence time was calculated using the formula T = Ks/2r Mya (Millions of years) where, Ks being the synonymous substitutions per site and r is considered 1.5 × 10-8 substitutions rate per site per year for dicot plants [44]. We reconstructed 24 conserved chromosomal blocks (labelled A-X) of B. rapa genome. Colour coding of the blocks was assigned based on their positions in a proposed ancestral karyotype (AK1-8) following the procedure stated by Cheng et al. [45] and Schranz et al. [46].

Microarray and microsynteny of PDI gene family
Two contrasting inbred Chinese cabbage (B. rapa ssp. pekinensis) lines, 'Chiifu' and 'Kenshin' in respect to cold tolerance and cold susceptible, were used for microarray expression analysis Kayum et al. [47]. A heat map was generated based on transcript abundance value of PDI genes using Cluster 3.0 (http://bonsai.hgc.jp/~mdehoon/ software/cluster/software.htm). The microsyntenic relationship of PDI genes among B. rapa, Brachypodium distachyon and A. thaliana were detected using Blast against whole genome of such crop species. Chromosomal positions of PDI genes were collected from respective databases and the relationship among the three crop species were plotted using circos software (http://circos.ca/) [48].

Plant growth conditions, treatments and sampling
Surface sterilize seeds of two Chinese cabbage inbred lines 'Chiifu' and 'Kenshin' were grown in incubation room maintaining 24°C temperature with 14/10 h (light/dark) condition at the Department of Horticulture, Sunchon National University, Korea. Abiotic stress treatments were imposed to four-week-old seedlings according to the methods described by Ahmed et al. [49], three biological replicates were maintained for abiotic stress treatments. Two contrasting Chinese cabbage inbred lines, cold tolerance 'Chiifu' and cold susceptible 'Kenshin' were used for abiotic stress cold and chilling injury stress experiment. Whereas, only 'Kenshin' was used to predict the BrPDI genes against other stress treatments (salt, drought, and ABA) due to unavailability of any contrasting genotypes for other abiotic stress treatments (salt, drought, and ABA) in our hand. By contrast, we assumed 'Kenshin' would be suitable for molecular characterizing of BrPDI genes for abiotic stress responsiveness rather than cold, because these type of genotypes are widely grown in tropics and sub-tropics. Fresh roots and leaves (third and fourth leaves) were excised from stress-treated plants over time points (0 h, 1 h, 4 h, 12 h, 24 h, and 48 h). Samples were immediately frozen in liquid nitrogen and stored at −80°C for RNA extraction. For chilling injury experiment, the Chinese cabbage inbred line 'Chiifu' and 'Kenshin' plants were grown in culture room under a 16 h light photoperiod maintaining 24°C. The three week old seedlings were transferred to incubator (TOGA clean system; model: TOGA UGSR01) maintained 0°C temperature and keep the seedling until clear remarkable cold injury symptoms appeared. Cold treated plants leaves were excised at 0 h, 24 h, 48 h and 72 h time points with the progression of cold injury in cold susceptible line 'Kenshin' (Fig. 9), thereby leaves were immediately frozen in liquid nitrogen and stored at −80°C for RNA extraction. Total five biological replications were maintained for chilling injury treatment. However, inbred line 'Chiifu' was used for organ-specific expression profiling of BrPDI genes and also infection with F. oxysporum f. sp. conglutinans for biotic stress as described by Ahmed et al. [50]. Fusarium-and mock-infected leaves (4th and 5th leaves) were collected at 0 h, 3 h, 6 h, 24 h, 3d, 6d, 9d, 12d after inoculation and immediately frozen in liquid nitrogen and stored at −80°C for RNA extraction.

RNA extraction and cDNA synthesis
Total RNA was extracted from the samples (roots and leaves) using an RNeasy mini kit (Qiagen, USA) following the manufacturer's instructions. DNA contamination was removed using RNase-free DNase (Promega, USA) following manufacturer's protocol. The extracted RNA was quantified by UV spectrophotometry at A260 using a NanoDropND-1000 and NanoDrop v3.7 software (Nano Drop Technologies, USA). Complementary DNA (cDNA) was synthesised from total RNA using a First-Strand cDNA synthesis kit (Invitrogen, Japan) following the manufacturer's instructions.

Qualitative and quantitative PCR expression analysis
Qualitative expression analysis was conducted using one-step Emerald Amp GT PCR Master Mix (Takara, Japan) by RT-PCR. Gene-specific primers for BrPDIs were used for RT-PCR and the BrActin gene of Chinese cabbage was used as an internal control. RT-PCR was performed using 1 μL of 50 ng template cDNA using master mix contained 10 pmol of each primer (forward and reverse), 9 μL sterile water and 8 μL Emerald mix in a total volume of 20 μL. The PCR conditions were initial denaturation 94°C for 5 min followed by 30 cycles of denaturation at 94°C for 30s, annealing at 58°C for 30s and extension at 72°C for 45 s, with a final extension at 72°C for 5 min. The PCR products were visualized on 1.2% agarose gel.
Quantitative Real-time PCR (qRT-PCR) was conducted using 10 μL reaction volume consisted of 5 μL 2× Quanti speed SYBR mix, 1 μL (10 pmol) each of forward (F) and reverse (R) gene-specific primers (Additional file 2: Table S3), 1 μL template cDNA (50 ng) and 2 μL ultra-pure water (ddH 2 O). The conditions for real-time PCR was initial denaturation at 95°C for 5 min, followed by 40 cycles of denaturation at 95°C for 10 s, annealing at 58°C for 10 s and extension at 72°C for 15 s. The qRT-PCR reactions were normalized using the Chinese cabbage Actin gene as a reference for all comparisons [51]. Fluorescence was measured at last step of each cycle, and three replicates were used for each sample. Amplification detection and data were processed using the Light cycler® 96 SW 1.1 software and the cq value was calculated using the 2 -ΔΔ C T method to determine the relative expression. The relative expression data was statistically analyzed (Tukey HSD test) and lettering was done using Minitab 17 software (https://www.minitab. com/products/minitab/).