Genome-Wide Identification and Analysis of Plant Cysteine Oxidase (PCO) Gene Family in Brassicaceae and Their Roles in Abiotic Stress Response

1. Introduction

Enzymes are the fundamental to all metabolic, catabolic and adaptive biological processes. They have evolved so that their functions meet biological needs, fine-tuning their activities to achieve optimal survival and replication fitness of the organism. Plant Cysteine Oxidases (PCOs) use molecular oxygen atoms to catalyze dioxygenation of Nt-Cys to Cys-sulfinic acid in ERF-VIIs (ETHYLENE RESPONSE FACTOR group VIIs, as well as other substrates VRN2 and ZPR2) [1,2] for subsequent Nt-arginylation [3,4]. Meanwhile, ERF-VIIs are known as the regulators of hypoxia-regulated transcriptional reprogramming to adapt the environmental change [5,6,7]. PCO activity is sensitive to physiologically relevant fluctuations in O₂ availability [8], thus these enzymes can act as plant O₂ sensors with a key role in regulating ERF-VII stability. In other words, as an enzyme, Plant Cysteine Oxidases (PCOs) are a direct link between environmental stimuli and molecular physiological outcomes [9].

In Arabidopsis, AtPCOs catalyze oxidation of cysteine to Cys-sulfinic acid at the N-termini of target proteins, a reaction that co-translational methionine cleavage exposes the N-terminal Cys for oxidation, leading to the degradation of the target protein via the Cys/Arg branch of the N-degron pathway [3,4,8]. Oxidized N-terminal Cys residues are substrates for arginyl transferase enzymes, with the arising arginylated N-termini recognized by ubiquitin ligases [10]. Ubiquitination is the signal for protein degradation by the proteasome. It is reported that PCOs isolated from Marchantia polymorpha (MpPCO) and Klebsormidium nitens (KnPCO) exhibit cysteine dioxygenase activity, indicating that PCO enzymes are conserved in early land plants and algae [11]. The MpPCO was incubated with the N-termini of MpERF-like (a 14-mer peptide representing the Cys-initiating N-terminus of MpERF-like, CRMNKRLGKGETGL), and MpPCO-catalyzed MpERF-like oxidation reached 89.3% after 1h. Meanwhile, homologs of arginyl-tRNA transferase (ATE) and E3 N-recognin, PROTEOLYSIS (PRT) 6 can be found in Marchantia polymorpha, suggesting the conserved way of PCO in catalyzing oxidation of cysteine to Cys-sulfinic acid at the N-termini of target proteins.

As PCOs regulate the target protein levels in N-degron pathway, it is critical to uncover PCOs structures to manipulate their enzyme activities for crop improvement [12]. The structures of AtPCO4_1, AtPCO4_2 (two different structures of AtPCO4 from independent crystallization conditions) and AtPCO5 are resolved to 1.82-, 1.24- and 1.91-Å resolutions [9]. There is a core double-stranded beta-helix (DSBH) supporting three histidine residues to coordinate the active site metal ions involved in catalysis. To recognize the active site, the Tyr₁₈₂-Ser₁₈₃-Ser₁₈₄-Glu₁₈₅-His₁₈₆-Asp₁₈₇-Arg₁₈₈-His₁₈₉-Cys₁₉₀ fragment is characterized by targeted mutagenesis. Furthermore, AtPCO4 (H164D) or AtPCO4 (D176N) fails to reconstitute the wild-type phenotype in 4pco (quadruple pco1/2/4/5 mutant plants) mutant, while AtPCO4 rescues the phenotype of 4pco mutant. Therefore, the sites (amino acid 160-190 position) of PCO4 play an important role in the enzyme activity [9,13,14].

The enzyme activities of PCOs are involved in plant development and stress responses. Quadruple pco1/2/4/5 mutant plants, but not triple pco1/2/4 mutants, in A. thaliana generated by T-DNA insertion have previously been shown to cause developmental defects, including retarded growth, leaf serration, and male and female sterility, along with constitutive expression of low-oxygen response genes [2,4,9]. Overexpression of AtPCO1 or AtPCO2 in Arabidopsis show smaller leaves and decreased biomass [4]. It suggests that appropriate PCOs level or enzyme activity is the fundamental to maintain plant growth and development. The PCOs catalyze the reaction of oxygen with the conserved N-terminal cysteine of ERF-VIIs to form cysteine sulfinic acid triggering degradation under normal condition. PCO activity decreases under hypoxia condition and then the stabilized ERF-VIIs active the expression of anaerobic genes in response to stress [3,15]. Previous studies have shown that aerobic or hypoxia related genes are down-regulated in either AtPCO1 or AtPCO2 over-expressed plants, whereas anaerobic genes are constitutively up-regulated in quadruple pco1/2/4/5 mutants [10]. Compared with wide type, plants over-expressing AtPCO1 or AtPCO2 show sensitive to prolonged submergence stress [4]. Meanwhile, overexpression of the ERF-VII RAP2.12 in Arabidopsis results in improved tolerance to submergence and up-regulation of genes associated with the hypoxic response [7,12]. Therefore, PCOs, as an O₂ sensor, play a critical role in stress response, especially hypoxia stress.

About 7,500 years ago, a natural heterotetraploid Brassica napus (AACC, 2n = 38) was formed by the hybridization of Brassica rapa (2n = 20, AA) and Brassica oleracea (2n = 18, CC) [16,17,18]. Brassica napus is an important oil crop, with the largest planting area as well as total oil production in China. With the changes of the global climate, the cultivation of rapeseed has been effected by various abiotic stress, such as flooding, waterlogging, drought and cold. Research on genetics, evolution and stress response of rapeseed is a good way to improve the quality and yield of rapeseed [19,20,21,22,23]. While PCOs are expected to play important and diverse roles in stress responses, detailed genome-wide analysis of PCO gene family has not been performed. In this study, we investigated the significant role of PCO genes in Brassicaceae plants. As a result, we identified 35 PCO genes of Brassica napus, Brassica rapa and Brassica oleracea and compared their gene structures, chromosomal locations, evolutionary relationships and expression patterns in different tissues and under different abiotic stress treatments. Therefore, the genome-wide analysis in this study provides a comprehensive understanding of PCOs in development and stress responses in Brassicaceae plants.

2. Results

2.1. Identification and Classification of PCO Genes in B. napus, B. oleracea and B. rapa

To identify PCO proteins in B. napus, B. oleracea and B. rapa, we performed a BLASTp search against the annotated proteins of B. napus, B. oleracea and B. rapa in BRAD (http://www.brassicadb.cn/#/BLAST/) using Arabidopsis AtPCO protein (AtPCO1, AtPCO2, AtPCO3, AtPCO4 and AtPCO5) sequences as queries. Sets of 20 (4 BnaPCO1, 2 BnaPCO2, 6 BnaPCO3, 4 BnaPCO4 and 4 BnaPCO5), 8 (2 BoPCO1, 1 BoPCO2, 1 BoPCO3, 2 BoPCO4 and 2 BoPCO5) and 7 (2 BraPCO1, 1 BraPCO2, 1 BraPCO3, 1 BraPCO4 and 2 BraPCO5) PCO proteins were identified in B. napus, B. oleracea and B. rapa, respectively (Table 1). All 35 PCOs in B. napus, B. oleracea and B. rapa encoded the amino acids ranging in length from 82 to 308 with the molecular weight from 8,650 to 34,180 Dalton. Isoelectric point (pI) of these amino acids ranged from 4.39 to 8.78 (Table 1).

To explore the classification and evolutionary characteristics of the PCO proteins, an unrooted phylogenetic tree based on the 40 protein sequences of B. napus (20), B. oleracea (8) and B. rapa (7) and Arabidopsis (5) PCO genes was constructed in MEGA 11 (Figure 1). According to the phylogenetic analysis, the PCO genes were divided into five groups: PCO1 (homologous to AT5G15120.1/AtPCO1), PCO2 (homologous to AT5G39890.1/AtPCO2), PCO3 (homologous to AT1G18490.1/AtPCO3), PCO4 (homologous to AT2G42670.1/AtPCO4) and PCO5 (homologous to AT3G58670.1/AtPCO5). In addition, according to the phylogenetic tree, PCOs can be clearly divided into three clades: PCO1s and PCO2s are clade A, PCO3s are clade B, and PCO4 and PCO5 are clade C, respectively.

2.2. Chromosomal Distribution and Duplication of BnaPCOs

As shown in Figure 2, the 20 BnaPCOs were scattered on 11 of the 19 B. napus chromosomes. The numbers of BnaPCOs on An-subgenome (8 genes) and Cn-subgenome (12 genes) showed a biased trend, with more genes on Cn-subgenome. B. napus (AnAnCnCn, 2n = 38) originated from the hybridization event between Brassica rapa (AnAn, n = 10) and Brassica oleracea (CnCn, n = 9) about 7500 years ago, and Brassicaceae species underwent a whole-genome triplication event [16]. Therefore, theoretically, B. rapa, B. oleracea, and B. napus genomes may have 15, 15, and 30 PCOs expanded from the 5 AtPCOs. In fact, only 7, 8 and 20 genes were identified in B. rapa, B. oleracea and B. napus, respectively (Table S1), indicating that some PCOs may have been lost during evolution.

The An and Cn subgenomes of B. napus were collinear with the corresponding diploid Ar and Co genomes, most of the An-Ar and Cn-Co homologous pairs showed similar chromosomal locations [16]. However, the C genome of B. napus had more homologous genes corresponding to the B. oleracea PCO3 gene (Bo5G025500), and the A genome had more homologous genes corresponding to the B. rapa PCO4 genes (BraA05G003220), while the homologous genes in the PCO1, PCO2 and PCO5 subgroups were relatively conserved. It was inferred that the tandem genome was mutated or duplicated in the homologous evolution process.

2.3. Gene Structures and Motif Analysis of PCOs in B. napus

In order to study the homology domain and conservation degree of the BnaPCOs, MEME [24] and TBtools [25] were used to predict and visualize their conserved domain and gene organization, respectively. 8 motifs were predicted by MEME (https://meme-suite.org/meme/tools/meme) (Figure S1). As shown in Figure 3, Motif 2 and 5 were the most conserved, with 38 of the 40 genes containing these two motifs. All of the PCO4 and PCO5 proteins contained Motif 8, whereas none of the other genes contained Motif 8. All PCO5 proteins contained eight motifs, of which the Motif 7 was located at the C terminal. However, all PCO2 proteins contained seven motifs in addition to the Motif 8, with the Motif 7 located at the N terminal. Two of the four short abnormal genes, BnaC05G0471900ZS and BnaC09G0521600ZS, lacked Motif 1-3-4-6-7-8, while BnaC04G0074400ZS and BnaC05G0283800ZS lacked Motif 1-2-4-5-6-8.

It is reported that some introns play an essential role in transcriptional regulation [26]. We also investigated the distribution of introns and exons to study the diversity of gene structure. Five exons and four introns were distributed on most of PCOs genes. And the intron phases of BnaPCOs were highly conserved in the same group, implicating the evolutionary similarity between these members.

2.4. Cis-Element Analysis of BnaPCOs

Cis-elements regulate the initiation and efficiency of gene transcription by binding to transcription factors [27]. We analyzed the cis-elements of 20 PCOs promoters with PlantCARE [28]. Most cis-elements of BnaPCOs promoters were associated with plant growth and development, hormone response and abiotic stress response (Table S2). Some cis-elements, such as, anaerobic induction, circadian control, defense and stress responsive, drought induction, light response, low temperature response and meristem expression, were the main members in regulating environmental stress. Among them, the most common elements were associated with light response and anaerobic induction, indicating that the growth and development of plants regulated by BnaPCOs was affected by light and oxygen. PCOs was known as a sensor of oxygen, which was consistent with the cis-elements result. However, the light response of PCOs needed further study.

In addition, there were more anaerobic-induction elements on PCO2s promoter region compared with other genes and there was no anaerobic-induction elements on the promoter of BnaC05G0471900ZS (belongs to BnaPCO3) (Figure 4). The cis-elements in the upstream promoter region of genes are closely related to the expression and function of downstream genes [29]. Therefore, according to the anaerobic-induction elements result, it indicated that the expression of PCO2s may be highly induced by hypoxia and it also implies the diversification of biological functions of PCOs genes in Brassica napus.

2.5. Expression Profiling of PCO Genes in Different Tissues

Based on Arabidopsis eFP Browser data (http://bar.utoronto.ca/efp/cgi-bin/efpWeb.cgi) and RNA-Seq data (B. rapa: GSE43245, B. oleateaca: GSE42891 and B. napus: BnIR database) (Table S3), the PCOs gene were expressed in different vegetative and reproductive organs of the four species at different developmental stages. In general, the expression pattern of PCO differed between groups (Figure 5). Almost all PCO were weakly expressed in pollen, indicating that PCO expression was down-regulated in sperm cells, presumably due to chromosomal structure or histone modifications. On the other hand, PCO5s were expressed at similar levels in all tissues, indicating that it appeared to be a housekeeping gene involved in plant growth and development.

There were six PCO3 genes in B. napus, three of which were expressed at high levels and the other three at low levels. The genes with higher expression levels were generated from B. oleateaca and B. rapa, and the expression pattern of these genes was also similar to their original genes. However, the genes with lower expression levels were extra mutated. This suggested that the previous derived gene may play a major role, while the new genes may be pseudogenes that have been amplified during evolution or may be silent under normal conditions and be as a backup for their homolog genes under special conditions.

2.6. Expression Profiling of PCO Genes under Abiotic Stress Treatment

To reveal the roles of PCOs in stress response in B. napus, the expression patterns upon various abiotic treatment were investigated (Figure 6, Table S4). In general, the expression of most PCO genes in leaves did not change significantly under various abiotic stress treatments (Figure 6a). In leaves, PCO3 gene expression was significantly induced under freezing stress, and PCO5 gene expression was increased under salt and osmotic stress as well. However, the response of PCO genes to stress was much stronger in roots than in leaves (Figure 6b). In roots, the expression of PCO1 and PCO2 were strongly induced by drought, and the expression of PCO1, PCO2, PCO3 and PCO4 were extremely down-regulated by freezing and cold stress. But PCO5 gene expression in roots did not change significantly under the abiotic stress treatments nonetheless, which was consistent with the hypothesis that PCO5 worked as a housekeeping gene.

Though there were similar motifs of drought-inducible and low-temperature responsive on PCOs promoter region, the gene expression pattern was various under drought, freezing or cold stress treatment. It implied that there were other regulators along with cis-elements to regulate PCO gene expression. It was interesting that PCO3 was up-regulated in leaves and down-regulated in roots after freezing treatment. In other words, PCO3 showed opposite response patterns in leaves and roots under freezing stress treatment. This suggested that PCO3 may play different roles in leaves and roots under freezing stress, and further studies were needed.

2.7. Expression Profiling of PCO Genes under Waterlogging Stress

Waterlogging removes air from soil leading to a blockage of gas exchange between the soil and the atmosphere [30,31]. Meanwhile, the diffusion rate of oxygen in water is only one tenth of that in air. As a result, oxygen availability in flooded soils is greatly limited, leading to suppression of root respiration. As mentioned before, PCO is an oxygen sensor in plant [8]. Therefore, it is speculated that PCO plays an important role in the hypoxic response induced by waterlogging stress. In order to elucidate the potential function of BnaPCO in response to waterlogging stress, the transcriptional patterns of BnaPCO in B. napus were studied after waterlogging stress. Transcriptome data (Figure 7, Table S5) showed that compared with CK, the expression of PCOs was significantly induced by waterlogging stress. Compared to other genes, the expression of PCO2 gene pairs (BnaA04G0111100ZS and BnaC04G0395100ZS) were strongly induced both in leaves and roots. While, the expression of PCO4 gene pairs (BnaA05G0034400ZS, BnaC03G0244400ZS and BnaC04G0037200ZS) were particularly strongly expressed in leaves. In addition, gene expression of PCO3 (BnaC05G0471900ZS and BnaC09G0521600ZS) were barely detectable in the transcriptome data, which was consistent with the previous result (Figure 5 and Figure 6). According to the expression results, the function of PCO was conserved under hypoxic stress. It suggested that BnaPCO was a vital component, connecting environmental stimulus with cellular and physiological response and BnaPCO could be a potential target for improving waterlogging stress tolerance.

3. Discussion

Plant cysteine oxidase family (PCO) is a set of plant O₂-sensing enzymes, which catalyze the O₂-dependent step [8]. In Arabidopsis, PCO has five members, PCO1, PCO2, PCO3, PCO4 and PCO5. In this study, we found 20, 8 and 7 PCO genes in B. napus, B. oleracea and B. rapa, respectively (Figure 1). In B. napus, the number of PCO genes in the An subgenome (8) was almost the same as that in their diploid ancestors B. rapa (7) (Table 1). This showed that all An subgenome PCO genes were preserved after the whole genome duplication event in B. napus. However, PCO3 genes in the Cn subgenome (BnaC04G0074400ZS, BnaC05G0155300ZS, BnaC05G0283800ZS, BnaC05G0471900ZS and BnaC09G0521600ZS) were mutated from the diploid ancestor B. oleracea (Bo5G025500). It suggested that Cn subgenome was much flexible during the evolution, compared with An subgenome [32]. PCO sequence alignment (Figure S2-S6) revealed that most PCOs were conserved in Brassicaceae, indicating that these duplicated PCO genes can still retain the function of their ancestors in these species.

The gene expression pattern of duplicated genes with similar functions may change during the formation of allopolyploids, which takes several typical patterns, including transgressive up−/downregulation, unequal parental contributions, and silencing [32,33]. Although PCO3 was expanded in B. napus, half of them were slightly expressed in different tissues or abiotic stress treatment. According to the results of PCO3 alignment (Figure S4), BnaA06G0127100ZS was derived from BraA06G014250, BnaC05G0155300ZS was derived from Bo5G025500, and the remaining four genes may be generated by mutation or amplification of DNA fragments during the evolution of polyploid. The expression patterns of BnaA06G0127100ZS and BnaC05G0155300ZS maintained their expression patterns in two diploid progenitors. In addition, the expression of BnaC04G0074400ZS and BnaC05G0283800ZS responded to abiotic stress (Figure 6). It indicated that the newly generated genes may contribute to phenotypic differences between allopolyploids and their parental species under abiotic stress conditions.

The binding between transcriptional factors (TFs) and cis-regulatory DNA sequences in the 5′-upstream of genes either activates or represses gene expression [34,35]. In general, genes with related expression patterns contain stress response elements in their promoter region and are likely to be regulated by stress related TFs [20,26,31,36,37]. It is reported that HSFB2b directly binds to GmC4H and GmCHS3 to regulate the gene expression in response to salt stress, since there are HSEs (Heat Shock Elements) in the promoter regions of the GmC4H and GmCHS3 [38]. It was showed that hormone-responsive elements and environmental stress-related elements were distributed on the BnaPCO promoters (Table S2). Combined with the expression data with abiotic stress treatment, BnaPCO expression was regulated by various stress responses, especially drought, freezing, cold and waterlogging stress (Figure 6 and 7). However, the expression level was different under different stress conditions. For instance, there were 7 anaerobic-induction cis-elements on the promoter region of BnaC04G0395100ZS, but the degree of waterlogging induced expression was much different between leaves and roots (Figure 7). As chromosome structure, histone modification, DNA methylation, transcriptional factors, cis-elements and other regulators work together to regulate gene expression [35], more studies are needed to explore the regulation mechanism of PCO expression under abiotic stress.

Oxygen homeostasis is critical for crop development, and hypoxia in plants is typically a consequence of reduced O₂ diffusion under conditions of waterlogging or submergence [19,39,40,41]. The quality and yield of rapeseed are seriously affected by waterlogging stress in China [23,41]. It is reported that the response to hypoxia in rice, Arabidopsis and barley is mediated by the group VII ETHYLENE RESPONSE FACTORs (ERF-VIIs) [1,5,6,7,15]. Moreover, the PCOs directly link O₂ availability to ERF-VII stability and anaerobic adaptation, leading to the suggestion that they act as plant O₂ sensors [3,4,8]. There were multiple anaerobic-induction cis-elements distributed on the BnaPCO promoters, and BnaPCO expression were induced by waterlogging stress after 6h treatment both in leaves and roots. Furthermore, H164 and D176 were in the AtPCO4 active site [9] and it was conserved in BnaPCO4 as shown in Figure S7. It implied that the function of BnaPCO in catalyzing ERF-VIIs could be conserved and targeting PCOs will be an effective way to improve the rapeseed tolerance to waterlogging stress by manipulating their O₂ sensitivity and/or substrate specificity.

4. Materials and Methods

4.1. Identification of the PCO Gene Family

The protein and nucleotide sequences of AtPCOs were obtained with TAIR (https://www.arabidopsis.org/). AtPCO proteins were used as query sequences to search for the PCO proteins of B. napus, B. rapa and B. oleracea using BLASTp (E-value < 1 × 10−5) in BARD (http://www.brassicadb.cn/#/BLAST/; protein databases were Brana ZS V2.0 pep, Brara Chiifu V3.5 pep and Braol JZS V2.0 pep, respectively). ExPASy (Gasteiger et al., 2003) was used to investigate the physical and chemical properties of these PCO proteins.

4.2. Phylogenetic Analysis and Chromosomal Locations

ClustalW was used to align the multiple sequences of all PCO proteins (from Arabidopsis, B. napus, B. rapa and B. oleracea), and MEGA 11 was used to build a phylogenetic tree using the neighbor-joining (NJ) phylogenetic technique with 1000 bootstrap replicates. The genome data and annotation information (gff3 files) were downloaded from public databases. TBtools version 1.116 was used to examine the gene chromosomal localization with the genome data and gff3 files [25].

4.3. Analysis of Gene Structure, Motif Composition and Cis-element

The motifs of BnaPCO proteins were predicted using the MEME v5.5.2 [24]. The number of motifs should not exceed 8. The distribution of motifs occurs zero or one time in each sequence. In order to investigate the structural characteristics of the BnaPCOs, the gff3 file was downloaded from the B. napus database (http://cbi.hzau.edu.cn/cgi-bin/rape/download_ext), which has the annotation information of the B. napus genome. The position information of introns and exons is obtained from the gff3 file. Meanwhile, the motif information was submitted to TBtools to graphically display gene structures and motif distributions.

To identify the cis-element of BnaPCOs, TBtools was used to obtain the 2,000 bp sequences in front of the genomic CDS. Then, the PlantCARE [28] was used to predict the cis- elements on these promoters. Thus, the number and types of different cis-acting elements in BnaPCOs were classified and visualized with TBtools.

4.4. Plant Materials and Treatments, Heat Map Analysis of the PCO Transcriptome Data

Zhongshuang 11 (the semi-winter cultivar of Brassica napus) seeds were germinated on filter paper, and the seedlings were then transplanted into pots with soil and nurtured in a growth chamber for four weeks (23 °C, 16 h light/8 h dark cycle). The waterlogging treatment was performed as previously described [22]. The pots of 4-week-old seedlings were placed in a 28cm x 14cm x 14cm container filled with water and the water level was maintained at about 2cm above the soil surface. Control plants (CK) remained well-watered throughout the experiment. Leaf and root samples of CK, 6h-treatment and 12h-treatment were collected for RNA isolation and RNA Sequencing was performed by the Illumina NovaSeq 6000 (Novogene). Hisat2 v2.0.5 was used to map the reads to reference genome and gene expression level was determined by FPKM calculation.

The expression patterns of BnaPCOs in different tissues and other abiotic stress treatment were obtained from the BnTIR (Brassica napus transcriptome information resource) database. The expression data of AtPCOs, BraPCOs and BoPCOs were obtained from Arabidopsis eFP Browser data, RNA-Seq data of B.rapa (GSE43245) and RNA-Seq data of B.oleareaca (GSE42891). All the expression data was standardized based on log2 scale, and clustered and visualized with TBtools.

5. Conclusions

In this study, we identified 20, 7 and 8 PCO (Plant Cysteine Oxidase) proteins in B. napus, B. rapa and B. oleracea, respectively, by exploring the important role of PCO genes in Brassicaceae plants. Collinearity analysis shows that the PCO gene family was relatively conserved in evolution of B. rapa, B. oleracea, Arabidopsis, and B. napus. However, there were mutations or duplications of PCO3 and PCO4 during homologous evolution process. The cis-elements that regulate plant growth and development and response to biotic and abiotic stresses were found in the BnaPCO promoters. In addition, we found that the Bna/Bra/BoPCO genes were expressed differently in different tissues at different developmental stages. Remarkably, some of BnaPCOs were significantly induced after waterlogging treatment, which was consistent with the cis-element analysis and previous studies. This study provides a foundation for further understanding the biology and stress response functions of PCO family genes in B. napus.