Genome-wide identification and transcriptomic data exploring of the cytochrome P450 family in Chinese cabbage (Brassica rapa L. ssp. pekinensis)

ABSTRACT Cytochrome P450 genes are involved in the catalysis of various reactions, including growth, development and secondary metabolite biosynthetic pathways. However, little is known about the characteristics and functions of the P450 gene family in Chinese cabbage (Brassica rapa L. ssp. pekinensis). In this study, we identified 258 non-redundant members of the cytochrome P450 gene superfamily in Brassica. rapa from Brassica database (http://www.brassicadb.cn). We also did transcriptomic data mining derived from abiotic/biotic stresses and developmental related RNA-seq in Chinese cabbage. These differentially expressed genes (DEGs) were involved in cadmium (Cd) detoxification, salinity tolerance, calcium signaling transduction and homeostasis, pathogen triggered immunity and leafy head formation in Chinese cabbage. These P450 genes are essential candidates for a better understanding of the mechanism of P450 genes, which help Chinese cabbage adapt to environmental changes.

Cytochrome P450 (P450s or CYP) belongs to a superfamily of ubiquitous heme-dependent enzymes in Bacteria, Archaea and Eukarya (Nelson et al. 2013). P450s are most well known in steroid/sterol biosynthesis and degradation, vitamin biosynthesis, xenobiotic detoxification, and metabolism of other drugs in plants and humans and play essential roles in the biosynthesis of vast arrays of secondary metabolites (Guengerich and Munro 2013). Plants have far more P450-encoding genes than mammals, and Arabidopsis thaliana has 249 and wheat has 1476 P450 genes, respectively (Guengerich 2019). P450 proteins' biological function in plants varies. These functions are involved in defence-related compounds, pigments biosynthesis, antioxidants and structural polymers etc. Besides, some P450s catalyze the metabolism of endogenous compounds, particularly those involved in signaling transduction. Therefore, P450s are essential bridge elements in mediating interactions between the organism and its chemical environment and regulating physiological processes.
Brassica rapa L. ssp. pekinensis belongs to the family of Brassicaceae, which is one of the most highly consumed vegetables throughout Eastern Asia. Besides, Brassica rapa is considered a model plant for evolutionary research on genome polyploidy (genome AA). Compared with the model plant Arabidopsis thaliana, Brassica species experienced an extra whole-genome triplication (WGT) event, which plays a pivotal role in the speciation and morphotype diversification of Brassica plants (Cheng et al. 2014). WGT results in a bulked gene pool that allows multicopy genes to evolve different functions or new functions, like subfunctionalization or neofunctionalization. According to evolutionary research, plant P450s belongs to two groups: A type and non-A type. P450s grouped to A type are specific to plants and have evolved since the divergence of plants from fungi and animals; while P450s grouped to non-A type are non-specific to plants but share sequence similarity to animals and fungal enzyme (Durst and Nelson 1995;Yu et al. 2017). The plant P450 families are grouped into nine different clans, in which the CYP71 clan was classified into Atype P450s, and the remaining eight clans, including CYP51, CYP710, CYP711, CYP72, CYP74, CYP85, CYP86 and CYP97 clans are grouped to non-A type (Nelson and Werck-Reichhart 2011). These clues contribute essential information to study the phylogeny of the CYP gene family.
In this study, we identified members of the cytochrome P450 gene superfamily in Brassica. rapa from Brassica database (http://www.brassicadb.cn) and compared with those in Arabidopsis. thaliana and rice-based on their phylogeny relationship. We also did a genome-wide bioinformatic analysis of P450 gene structure and motifs, genomic distribution, possible SSR marker loci and chromosome collinearity analysis. Besides, we did RNA-seq data mining from different stress treatment to unveil the transcriptome   dynamics of Brassica rapa P450 genes and identified the candidate P450 genes that are conserved, specific, or show correlated expression with these abiotic stresses. We hope the detailed information from this research could further investigate the P450 genes and their specific function in Chinese cabbage.

Identification of the Chinese cabbage P450s and evolutionary relationship
Based on the Brassica. rapa line Chiifu's genomic and annotated genes information, after filtering the published Bra040154 is the shortest (162 bp), and Bra028715 is the longest (3882 bp). The proteins encoded by the BrP450s were predicted to show significant differences in their sizes and physicochemical properties (Table 1). The predicted protein lengths varied from 53 to 1293 amino acids. Besides, P450 protein sequences had significant variations in predicted isoelectric point (pI) values (ranging from 4.9 to 11.38) and molecular weights (ranging from 5.99 to 61.45 kDa). Two hundred fifty-eight BrP450 genes, two hundred fifty-one, could be unevenly distributed on ten chromosomes and seven located on six Scaffolds (Figure 1). Chromosome one carries twentyfour BrP450s, chromosome two has twenty-nine, whereas chromosomes three, four, five, six, seven, eight, nine and ten each contain thirty-six, eighteen, twenty, thirty-five, thirteen, nineteen, thirty-nine and fourteen, respectively.
Gene structure and motif divergence of BrP450 genes Gene structure and motifs often diverge during the evolution of multi-gene families. To explore the phylogenetic relationship and pattern of the BrP450 gene structures, we analyzed the intron-exon and motif characteristics of BrP450s by aligning the CDSs to the genomic sequence. The phylogenetic relationship of BrP450 was shown in (Figure 2A). The identified BrP450 genes have a minimum zero and maximum of 11 introns ( Figure 2B). Out of the 258 identified BrP450 genes, 50 (19.38%) genes have no intron, 93 (36.05%) genes have one intron, 39 (15.12%) genes have two introns, 19 (7.36%) genes have three introns, 17 (6.59%) genes have four introns, and 40 (15.50%) genes have five or more than five introns. Additionally, the number of introns in the CDSs of P450 genes in the same subfamily was different.
Simple sequence repeat (SSR) markers are a type of codominant marker and have many advantages like multi-allelic variation, abundance throughout the genome, a high level of polymorphism, etc. So it is suitable to develop molecular markers that associated with genetic mapping and fingerprint assistant plant breeding (Powell et al. 1996). In this study, 85 SSR markers were detected in the 258 identified BrP450s using the online SSR identification tool SSRIT, which all belong to the di-motif type (Table 2). Besides, 69 out of 85 detected with SSR marker P450 genes have only one SSR marker, and 14 out of 85 P450 genes have two SSR markers, and 2 out of 85 P450 genes have three SSR markers. Of the detected SSRs loci in the BrP450 genes, 76 were detected in exons, and 27 were detected in introns.

Conserved domains and motifs in the predicted BrP450 proteins
The MEME web server was used to analyze motif distribution and to conduct domain predictions. A total of 15 motifs were identified in the BrP450s ( Figure 3A, B and C). The motifs for BrP450s belonging to the same cluster were conserved. To analyze the motifs' contribution to BrP450s gene family, we calculated the frequency of these motifs.   To gain insight into the evolution of P450s in Chinese cabbage, we constructed a Neighbor-Joining phylogenetic tree using the P450 proteins from representation plant species, including the rice P450s (OsP450) and Arabidopsis P450s (AtP450) (Figure 4). The 258 BrP450 proteins were classified into A-type and non-A type P450 gene families. The non-A type gene family contain eight major clusters, designated as CYP72 clan, CYP97 clan, CYP711 clan, CYP710 clan, CYP51 clan, CYP74 clan, CYP85 clan and CYP86 clan. The A-type P450 gene family were clustered into the CYP71 clan. The phylogenetic tree suggested that the BrP450 proteins were closest to the AtP450 proteins and more distance from the OsP450 proteins.

Differential expression profiles of BrP450 genes
To indicate the potential functions of P450 genes in Chinese cabbage, we analyzed stresses-treated RNA-Seq data from expression profiles of different biotic and abiotic stress. The BrP450s exhibited a diverse expression pattern. Among the 194 Cd triggered DEGs, 11, 18 and 12 genes were up-regulated, and 15, 8 and 9 genes were down-regulated in 25, 50 and 100 mM Cd treatment, respectively  (Figure 7a; Table S1). We also conducted a cluster analysis for these DEGs and obtained robust cluster results for the up-and down-regulated genes ( Figure 7a). Besides, the Veen diagram indicated there are 7 P450 DEGs were identified co-expression at each Cd concentration (Figure 7b). Among these genes, 6 DEGs were up-regulated, and one DEG was down-regulated. We also assigned the commonly expressed P450 genes to GO terms by blast2go ( Figure S1). Most DEGs enriched subcategories belonged to the categories of biological process (2 subcategories), cellular component (4 subcategories) and molecular function (6 subcategories). Based on this study, these P450 genes were inferred to be involved in the stress response to Cd in Chinese cabbage, and they may have important regulatory roles in the Cd detoxification for Chinese cabbage.
To explore the salinity response P450 genes, we found two DEGs between the salinity treatment group and the CK group. Interesting, these two genes were all down-regulated in the salinity treatment group compared with the CK group ( Figure 8; Table S2). The expressed P450 genes to GO terms were enriched by blast2go ( Figure S2). This result indicates that these genes may involve negative regulation in salinity tolerance.
To identify the calcium deficiency response P450 genes, which could be the possible reason for Chinese cabbage tip burn symptom, 7 DEGs were found between the calcium deficiency treatment group and the CK group. Among these genes, five genes were up-regulated, and two genes were minor down-regulated in calcium deficiency treatment compared with CK ( Figure 9; Table S3). The expressed P450 genes under calcium deficiency treatment to GO terms were enriched by blast2go ( Figure S3). This result showed that these P450 genes might participate in the plant calcium signaling system to maintain the Chinese cabbage cytosolic calcium homeostasis.
This study also analyzed the transcriptome data to indicate biotic stress response P450 genes caused by downy mildew. There are 9 DEGs found between the 13-318-R-CK (resistant variety) group and 13-361-S-CK (susceptible variety) group. Among these genes, 3 DEGs were up-regulated in the 13-361-S-CK group compared with the 13-318-R-CK group, and 5 DEGs were down-regulated in the 13-361-S-CK group compared with 13-318-R-CK group ( Figure  10a; Table S5). There are 12 DEGs found between the 13-318-R-CK group and 13-318-R-T (treated by downy mildew) group. Among these genes, 6 DEGs were up-regulated in the 13-318-R-T group compared with the 13-318-R-CK group, and only one DEG was down-regulated in the 13-318-R-T group compared with the 13-318-R-CK group (Figure 10b; Table S5). Besides, there are 10 DEGs found between the 13-318-R-T group and 13-361-S-T group. Among these genes, 4 DEGs were up-regulated in the 13-361-S-T group compared with the 13-318-R-T group, and one DEG was down-regulated in the 13-361-S-T group compared with the 13-318-R-T group (Figure 10c; Table S5). There are also 16 DEGs were found between the 13-361-S-CK group and 13-361-S-T group. Among these genes, 11 DEGs were up-regulated in the 13-361-S-T group compared with the 13-361-S-CK group, and there are no DEGs significantly (absolute log 2 FC value ≥1.5) down-regulated in the 13-361-S-T group compared with the 13-361-S-CK group (Figure 10d; Table S4). We also generated a Veen diagram to indicate the co-expression genes between each group ( Figure 10e). This indicates that these P450 genes may play an essential role in Chinese cabbage pathogen triggered immunity (PTI). To identify the correlation between P450 genes and Chinese cabbage leaf head development, we found 18 DEGs between the rosette leave (RL) stage and folding leave (RL) stage. Among these genes, five genes were up-regulated, and three genes were down-regulated in the folding leave stage compared with the rosette leave stage (Figure 11; Table S5). This result indicates that these P450 genes may involve in the Chinese cabbage leaf head development.
This research analyzed five kinds of abiotic and biotic stresses induced P450 genes in Chinese cabbage. We made a Veen diagram to identify the common expressed DEGs among the five stresses ( Figure S6; Table S6). Unfortunately, there are no common DEGs expressed among all of the five stresses.

Expression patterns of the BrP450 genes
BrP450 genes are essential in Chinese cabbage response to environmental abiotic and biotic stresses. To verify the expression pattern of the BrP450 DEGs derived from the RNA-seq, we used real-time quantitative PCR (RT-qPCR) for validation. The results indicated that the BraA01g001490, BraA06g008660, BraA08g031540 and BraA10g030360 had higher expression level on 50 mM Cd than the rest of Cd treatment, whereas BraA03g038810 had a higher expression level on 100 mM Cd than the rest of Cd treatment. Besides, BraA02g012430 and BraA03g031300 showed lower expression level in all Cd treatment than in control. (Figure 12).
The BrP450 response to salinity stress was also investigated. Compared with the control (CK), the transcript level of both BraA04g029510 and BraA09g053870 increased to about two-fold higher at 3 hours post-treatment. However, the transcript level on other time points of both genes showed lower than control (Figure 13).
The qPCR results indicated the BraA01g000830, BraA02g039080, BraA03g048740, BraA04g020810, BraA05g012440, and BraA10g003760 had lower expression level on calcium deficiency treatment than the control. Only BraA09g021180 had a higher expression level on 1/8 Ca treatment than the control (Figure 14).
The qPCR results showed different expression pattern under downy mildew infection. On the 13-13-3 infected plants, the BraA02g026450 and BraA05g010620 had higher expression level compared with control. While on 15-14-3 infected plants, the BraA05g005260 had a higher expression level compared with control. On the 15-44-3 infected plants, all of the seven BrP450 genes showed a higher expression   level than control. Moreover, on the 15-193-3 infected plants, the BraA02g026450, BraA05g005260, BraA05g010620, BraA06g037550, and BraA09g007830 had increased expression level compared with control. Besides, the BraA05g010620 and BraA06g037550 showed undetectable in 15-14-3 and 13-13-3 infected plants, respectively ( Figure 15).
During the Chinese cabbage leafy head formation process, the BrP450 also showed diverse expression patterns in different cultivars. In variety ST, six BrP450 genes except BraA09g003910 showed increased expression level in the folding stage, compared with the seedling stage. While in variety 495-4, only BraA02g041800 showed increased expression level, compared with the seedling stage. These results indicate that the BrP450 genes might have different expression pattern in different cultivars ( Figure 16).

Discussion
Cytochrome P450 genes are involved in the catalysis of various reactions, including growth, development and secondary metabolite biosynthetic pathways (Vasav and Barvkar 2019). As a critical supergene family, cytochrome P450s have been identified in various plants at the genome-wide level. However, little is known about these P450 genes' response to abiotic and biotic stress and how these P450 genes participate in leaf head development in Chinese cabbage. Besides, the second generation sequencing-based RNA-Seq technology is a powerful tool for detecting the whole-genome gene expression. Thus, a comprehensive analysis of the RNA-Seq data with specific environmental clues can help us to understand the mechanisms of P450 genes in Chinese cabbage. Here, we analyzed RNA-Seq data related to heavy metal (Cd) stress, salinity stress, downy mildew stress, tipburn and leaf head development. In this study, we identified 258 non-redundant out of 446 P450 genes from Chinese cabbage. Phylogenetic analysis of the P450s in different species indicated that the P450s in Chinese cabbage had 1.4-fold as many P450 members as Arabidopsis (Yu et al. 2017). Besides, each Chinese cabbage subcluster was closer to Arabidopsis than to the rice, which was consistent with the known evolutionary relationships of this species (Durst and Nelson 1995;Yu et al. 2017).
The Brassica lineage underwent WGT and generated triplicated genomic regions in Brassica rapa after the Brassica lineage split from a common ancestor with Arabidopsis thaliana.
Analysis of gene expression patterns can be used to predict the molecular functions of genes involved indifferent processes. Most cytochrome P450 genes induced by both abiotic and biotic stresses contained MYB and MYC's recognition sites, ACGT-core sequence, TGA-box and W-box for WRKY transcription factors in their promoters (Rudolf et al. 2017). These cis-acting elements are known to participate in the regulation of plant defence. The response of each gene to multiple stresses is strictly regulated. The plant P450 has two primary pathways: biosynthetic pathways, and another is involved in detoxification pathways.
The eighty-five SSR markers developed from this research might contribute to marker-assisstant plant breeding. SSR can be identified from either genomic DNA or cDNA sequences. BrP450s family genes abundance in the Brassica rapa genome, reproducibility and high polymorphic nature make BrP450s SSR markers more attractive and reliable. SSR markers derived from genomic libraries amplified more fragments and showed more polymorphism within crops. Besides, the BrP450 simple sequence repeats (SSRs) are repeat based and co-dominant markers. Di, tri and tetranucleotide repeats are extensively distributed all over the Chinese cabbage genomes. The microsatellite is a more powerful marker system due to the high level of allelic   variation. These are widely used due to DNA usage, low cost for manual assays and high throughput genotyping. SSRs have been widely used in developing genomic linkage maps, QTLs mapping and other genetic analysis of germplasm. Many plant breeders considered that SSRs are a more successful marker system for genetic diversity, identification of hybrids, classification and DNA fingerprinting of germplasm. SSRs have been extensively used for genetic analysis of many horticultural crops, including pear (Pyrus communis), grapes (Vitis vinifera), cucumber (Cucumis sativus), sugar beet (Beta vulgaris) and brassica species (Parthiban et al. 2018;Taheri et al. 2018;Wang et al. 2018a;Ahmad et al. 2020;Younis et al. 2020).
The plant detoxification mechanisms of heavy metals is a sophisticated regulation. Similar to other stresses, heavy metal toxicity induces typically oxidative stress and modulates the expression of various stress-related genes in plants (Shri et al. 2009;Preeti et al. 2012). Efflux of heavy metals from the plant roots could be an impending approach to reduce heavy metal build-up in plants. In this study, 4 out of 6 DEGs were enriched to the intrinsic component of the membrane (GO: 0031224), suggesting these BrP450 genes might involve signal recognition and transduction as reported in the previous research (Dubey et al. 2018). After chelation and sequestration of the heavy metal into subcellular compartments, the organic acid and amino acids will play a pivotal role in heavy metal detoxification. When the heavy metal ions were transported through the xylem and sequestrating ions in the vacuole, organic acids such as malate, citrate, and Besides, by metal chelation, amino acids and their derivatives also have an essential role in plant resistance to toxic levels of metal ions (Anjum et al. 2015). In our research, the DEG (BraA03g031300) derived from Cd stress was mapped to a predicted KEGG pathway (ko: 00380), which belonged to tryptophan metabolism. Plants possess various mechanisms through which they provide resistance against heavy metals toxicity, and antioxidant defence system that comprises antioxidant enzymes also have a significant contribution in detoxification (Yang and Chu 2011;Das and Roychoudhury 2014;Emamverdian et al. 2015). Among the DEGs, We found both BraA03g031300 and BraA03g038810 participated in the critical enzyme (abscisic acid 8 ′ -hydroxylase, EC 1. 14. 14. 137) biosynthesis in the oxidative catabolism of abscisic acid, which is considerable evidence to Chinese cabbage heavy metal stress response and detoxification. Transcriptomic analysis of Arabidopsis thaliana exposed to Hg reported that Hginduced genes encode proteins involved in P450-mediated biosynthesis of secondary metabolites, cell wall metabolism, and chlorophyl synthesis (Heidenreich 2001). In our study, BrP450 genes also participated steroid hormone biosynthesis (ko00140), retinol metabolism (ko00830), xenobiotics metabolism (ko00980), which might indicate the heavy metal response mechanism in Cruciferae species. Thus, data mining of these heavy metal-induced BrP450 genes could better understand the mechanism of tolerance and detoxification of Cd in Chinese cabbage.
Saline stress in the short term is displayed first in the roots by suffering osmotic stress associated with the accumulation of phytotoxic ions. In the long term, salinity caused ion toxicity to result in a cytosolic nutrient imbalance. The previous transcriptome analysis of short-period acclimation (for 24 h) to salt stress (200 mM NaCl) in the algae Chlamydomonas reinhardtii found that the regulation of phosphatidic acid acetate levels and lipid homeostasis had a crucial role in increasing tolerance to salt stress and use as alternative energy for solving the impairment of photosynthesis as Figure 11. Heatmap of P450 DEGs related to Chinese cabbage leafy head formation stages. well as the enhancement of glycolysis metabolism (Wang et al. 2018b). Another report shows that CrebZIPs TFs may play essential roles in mediating photosynthesis by reduced chlorophyl content and Fv/Fm ratio, while increased oil contents, carotenoid, and NPQ could be interpreted as adaptive mechanisms to salt stress (Ji et al. 2018). Research on Grapevine (Vitis vinifera) demonstrated, in response to salt stress, many genes were involved in various defence-related biological pathways, including ROS burst, heat shock proteins (HSPs), ion transport, pathogenesis-related proteins (PRs), and hormone signaling (Guan et al. 2018). In our study, the two BrP450 DEGs were enriched in the oxidationreduction process (GO: 0055114), and the KEGG enriched in glucosinolate biosynthesis (ko00966), 2-oxocarboxylic acid metabolism (ko01210) and biosynthesis of secondary metabolites (ko01110), which are classical pathways to saltstress in the previous report (Nazar et al. 2011;Aghajanzadeh et al. 2018). These results indicate the P450s also involved in Chinese cabbage salinity stress tolerance and resistance through a negative regulation loop.
The Chinese cabbage tipburn is thought to be a physiogene problem primary due to inadequate calcium distribution, and when it occurs, the inside of the leaf head became brown and then spread to the whole plant. Su et al. found that the cytoplasmic Ca 2+ fluctuation-induced downstream signaling events and SA signaling or other biological events are involved in the plant defence response to tipburn in Chinese cabbage (Su et al. 2016). Besides, some QTLs also found to enhance the Chinese cabbage resistance to tipburn (Su et al. 2019). However, the mechanisms of how Chinese cabbage response to this disease is still poorly understood. Our study found BraA05g012440 involved in the carotenoid biosynthesis pathway (CBP) (K04123). The carotenoid biosynthesis pathway is essential for producing  photosynthetic and protective pigments, plant hormones, and visual/olfactory attractants for animal pollinators and seed dispersers (Stanley and Yuan 2019). This result indicates that Chinese cabbage tipburn caused by calcium deficiency usually evokes a series of poor calcium distribution symptoms, including stress-related hormone pathways, like the ABA pathway and carotenoid biosynthesis pathway. Moreover, P450 genes might participate in these signal pathways in Chinese cabbage adaptation to a calcium-deficient environment.
Downy mildew is one of the most notorious diseases of Chinese cabbage caused by Hyaloperonospora brassicae. Genome-wide gene expression profiles in response to downy mildew in Chinese cabbage shows the differentially expressed genes identified between susceptible and resistant varieties are mainly involved in the plant-pathogen interaction pathway, plant hormone signal transduction pathway, and biosynthesis of secondary metabolites pathway . In this study, we found the P450 DEGs were involved in the oxidation-reduction process (GO: 0055114), which belong to the biological process category (GO: 0008150), and oxidoreductase activity (GO: 0016705), iron ion binding (GO: 0005506), heme binding (GO: 0020037), which enriched to molecular function category (GO: 0003674). This result suggests downy mildew activated P450 DEGs might be involved in pathogen triggered immunity pathway.
The leafy head formation of Chinese cabbage is thought to be determined by genetic factors and environmental clues. Genome-wide transcriptome analysis between Chinese cabbage and non-heading pak-choi revealed molecular pathways involved in leafy head formation (Sun et al. 2019). The previous research shows genes associated with photosynthesis-related processes such as chloroplast RNA binding in the RNA category, germin-like proteins (GLP3 and GLP1) in the stress category, isoprenoids carotenoids phytoene synthase (PSY), flavonoids dihydroflavonols (ATCHIL), and chloroplast lipoxygenase (LOX2) in hormone and secondary metabolism categories (Su et al. 2019). Recently released research reported, Chinese cabbage, BcpLH (Brassica rapa ssp. pekinensis LEAFY HEADS), a close homolog of HYL1 (HYPONASTIC LEAVES 1, Arabidopsis), is differentially expressed in juvenile leaves, which are flat, and in adult leaves, which display extreme incurvature, by affecting a specific subset of miRNAs and coordinates the direction, extent, and timing of leaf curvature during head formation period (Ren et al. 2020). In this study, we found the BrP450 DEGs primarily involved two GO categories, they are oxidation-reduction process (GO: 0055114) in the biological process category, and oxidoreductase activity (GO: 0016705), iron ion binding (GO: 0005506), and heme-binding (GO: 0020037) in molecular function category. Among the DGEs, we found Bra03623 involved in brassinosteroid-6-oxidase 1 biosynthesis pathway (K09590). Brassinosteroids (BRs) are steroid hormones that are essential for plant growth and development. These hormones control the division, elongation and differentiation of various cell types throughout the entire plant life cycle (Planas-Riverola et al. 2019). Abscisic acid (ABA) is involved in several critical processes in normal growth and development and adaptive responses to environmental stresses. This study, we found that Bra021965 participate in the abscisic acid 8 ′ -hydroxylase biosynthesis pathway (K09843). There is a possibility that the ABA spatiotemporal accumulation and distribution during the leafy head formation stage, which triggered a group of the downstream signal cascade.

Conclusion
In this research, we identified 258 non-redundant P450 genes in Chinese cabbage. The phylogenetic analysis suggested these BrP450s were more closely with AtP450s than with OsP450s. We also analyzed abiotic/biotic and developmental related P450 DEGs' expression pattern from transcriptomic data. These P450 DEGs were significantly enriched in Cd detoxification, salinity tolerance, calcium signaling transduction and homeostasis, pathogen triggered immunity and leafy head formation pathways in Chinese cabbage. This information will help better understand the P450 superfamily function in Chinese cabbage and provide potential molecular genetic targets for engineering stresses resistance crop in agricultural activity.

Identification and analysis of putative P450 genes in Chinese cabbage
Chinese cabbage genome annotation was downloaded from the Brassica database (BRAD, http://brassicadb.org/brad/), including the complete Brassica A genome sequence from B.rapa (Chiifu-401). The P450s nucleotide sequences were aligned using SnapGene5.0.5 Software (from GSL Biotech; available at snapgene.com). Schematic gene structure diagrams were drawn using Exon/intron organization of the P450 genes were determined using the Gene Structure Display Server (http://gsds.cbi.pku.edu.cn/) with GFF file. The motif of P450 protein sequences was analyzed with MEME online tools (http://meme-suite.org/tools/meme). For phylogenetic analysis, protein sequences of P450 from Chinese cabbage, Arabidopsis and rice were aligned with MUSCLE, a multiple sequence alignment tool from SnapGene5.0.5 Software (from GSL Biotech; available at snapgene.com) with default settings. The phylogenetic tree was constructed by the Neighbour-Joining (NJ) method and bootstrap test using MEGA7 Software with the following parameter settings: pairwise deletion mode, Poisson correction, and bootstrapping (1000 replicates). The molecular masses of putative P450 proteins were calculated using the Compute pI/Mw tool of ExPaSy (http://web.expasy.org/compute_pi/). The number of amino acids, molecular weight (MW), and theoretical isoelectric point (pI) were computed using the Prot-Param tool (https://web.expasy.org/protparam/). The SSR marker was detected using the SSRIT software (https://web. expasy.org/protparam/) with the parameters adjusted to identify perfect di-, tri, tetra-, penta-, and hexa-nucleotide motifs with a minimum of 6, 5, 5, 4 and 4 repeats, respectively. Based on the general feature format (GFF) file of the P450 genes in Brassica rapa L. ssp. Pekinensis. The physical map of P450 locate in chromosomes was generated by Map-Chart developed by Wageningen University and the research center (https://www.wur.nl/en/show/Mapchart.htm). The GO analysis of P450 proteins was carried out using the GO term analysis tool in Gramene (http://www.geneontology. org/) developed by the GO consortium. The DEGs and Veen diagram's cluster analysis was carried out by the TBtools Software (https://github.com/CJ-Chen/TBtools). Moreover, the commonly expressed P450 genes were assigned to GO terms by Blast2go Software. Conserved motifs in the full-length amino acid sequences of P450 proteins from Chinese cabbage were identified by MEME (Bailey and Elkan 1995). The Arabidopsis and rice P450 protein sequences were downloaded from the Arabidopsis Information Resource (http://www.arabidopsis.org) and the Institute for Genomic Research Rice Genome Annotation Project (TIGR: http://www.tigr.org/), respectively. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed to investigate the biological system's highlevel functions and utilities. The cluster Profiler R package was used to test the statistical enrichment of DEGs in KEGG pathways.

Plant materials growth condition
Chinese cabbage (Brassica rapa L. ssp. Pekinensis inbred line) plants were grown in the greenhouse at 20 ± 2°C with a photoperiod of 16 h light and 8 h dark.

RNA-seq data analysis of cytochrome P450s
To study gene expression profiles of the Chinese cabbage P450 gene family, we used RNA-seq data from previous studies. RNA-seq data from four different Cd concentration treatment (unpublished data), including 25, 50, 100 and 0 mM Cd, is a control in this experiment. The RNA-seq data of salinity stress using 200 mM sodium chloride treatment and 0 mM sodium chloride treatment is a control in the experiment (Qiu et al. 2017). The RNA-seq data of Chinese cabbage physiological caused tipburn using standard ½ MS medium (220 mg/L CaCl 2 •2H 2 O) and 0 mM Ca 2+ ½ MS medium, which is a treatment concentration (unpublished data). The RNA-seq data of downy mildew stress using Hyalo-peronospora brassicae to infect Chinese cabbage, and use 316S (susceptible line) and 318R (resistance line) as infected samples . The RNA-seq data of the leafy head formation experiment use rosette leave stage and folding leave stage as comparison samples (Wang et al. 2012). For the expression values of Brassica rapa L. ssp. Pekinensis cytochrome P450 genes, transcript abundance, were calculated by fragments per kilobase of exon model per million mapped reads (FPKM) from CuffLinks v2.1.0 the FPKM values were log2 transformed. A heat map was generated using MeV software after original the FPKM values were log2-transformed and centered. In this study, we use log2FC as a threshold, and genes expression log2FC more than 1.5 means up-regulated, and lower than minus 1.5 means down-regulated.

RNA extraction and real-time quantitative PCR
Total RNA was extracted from each sample using Trizol reagent (Invitrogen, Carlsbad, CA, USA) and treated with Rnase-free Dnase I (TaKaRa, Dalian, China) 45 min according to the manufacturer's protocol. First-strand cDNA was synthesized from 1 µg of total RNA using a PrimeScript 1st Strand cDNA Synthesis Kit (TaKaRa). RT-qPCR was carried out using SYBR Green Master Mix (TaKaRa) by IQ5 Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). The gene-specific primers designed for each gene are listed in Supplementary Data. The actin gene was used as a constitutive expression control in the RT-qPCR experiments. The PCR cycling conditions comprised an initial polymerase activation step of 95°C for 1 min, followed by 40 cycles of 95°C for 10 s and 60°C for 30 s. After each PCR run, a dissociation curve was designed to confirm the product's specificity and avoid primer dimers' production. The relative amounts of the amplification products were calculated by the comparative 2 −ΔΔCT method.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Notes on contributors
Shu Zhang is a PhD student at College of Life Sciences, Zhejiang University.
Qian-Rong Wu is a master student at School of Education, Huazhong University of Science and Technology.
Hui-Min Zhang is a master student at College of Life Sciences, Shandong Normal University.
Zhen-Ming Pei is an associate Professor and biologist in Duke University.