Genome-wide evolution and expression analysis of the MYB-CC gene family in Brassica spp.

Plant materials and growth conditions

B. napus (accession Zhongshuang11) was kindly provided by the Oilcrops Research Institute, Chinese Academy of Agricultural Sciences. The seeds were disinfected in 70% ethanol for 1 min, washed with ddH₂O three to four times, placed on a Petri dish lined with two layers of filter paper, and then germinated at 24 °C in growth chamber. After 3 days, the germinated seedlings were transferred to MS medium containing 0.3% agar and 3% sucrose at pH 5.7. The seedlings grew at 24 °C under a long-day light cycle (16 h light/8 h dark) and a relative humidity of 65%. Then the two-week-old plants were cultured in 1 mM and 1 μM Pi liquid medium for 24 h, respectively. Rapeseed seedlings’ shoots and roots were collected and frozen in liquid nitrogen. The RNA-seq analysis was carried out by Biomarker (Beijing). A total of 16 RNA samples were sequenced using the Illumina HiSeq 2,000 platform (Illumina, USA). The transcript abundance (FPKM value) was calculated based on the length of the gene and the reads mapped gene. The heatmap was constructed using TBtools (Chen et al., 2020). The raw RNA-seq data were deposited in the NCBI SRA database (https://ncbi.nlm.nih.gov/sra/) under the accession number PRJNA739537 and the assembled transcriptome data were deposited in the NCBI GEO database (https://www.ncbi.nlm.nih.gov/geo/)under the accession number GSE192452.

Identification of MYB-CC gene family in B. napus, B. rapa, and B. oleracea

The genes and proteins sequences in B. napus, B. rapa, and B. oleracea were downloaded from two databases (Genoscope: https://www.genoscope.cns.fr/brassicanapus/ and BRAD: http://brassicadb.cn/). The MYB-CC members were selected from B. napus, B. rapa, and B. oleracea by considering the E value < e⁻¹⁰ when compared with the 15 known MYB-CC protein sequences from A. thaliana (TAIR: http://www.arabidopsis.org/). This was achieved using BLAST on a computer (Rubio et al., 2001). The candidate protein sequences were further filtered using PFAM (http://pfam.xfam.org/) (Finn et al., 2014) and InterProScan (http://www.ebi.ac.uk/interpro) (Mitchell et al., 2015) according to the MYB and CC domains. The molecular weight and isoelectric point of MYB-CC members were calculated with ExPASy (Gasteiger et al., 2003).

Multiple sequence alignment and phylogenetic analysis of MYB-CC family

We constructed multiple sequence alignments of BnaMYB-CC, BraMYB-CC, BolMYB-CC, and AtMYB-CC using Clustal X (Thompson et al., 1997). The unrooted phylogenetic trees were constructed using the Neighbor-Joining (NJ) method. This was done using MEGA 6.0 software with the JTT model and pairwise gap deletion option (Tamura et al., 2013). The bootstrap analysis was conducted with 1,000 iterations.

Chromosomal location and gene duplication

The position of the MYB-CC genes was investigated according to the B. napus v4.1.chromosomal database (http://www.genoscope.cns.fr) and the Brassica database (BRAD, http://brassicadb.org). MapInspect software was used to locate the MYB-CC genes. The gene duplication was defined as: a coverage region between two sequences of more than 70% with more than 85% similarity in the coverage region (Zhou et al., 2004; Yang et al., 2008). Additionally, the ka (non-synonymous mutation rate) and ks (synonymous mutation rate) values of the repeating gene pairs were calculated using DnaSp software.

Gene structure and conserved motif analysis

GSDS (http://gsds.gao-lab.org/) was used to determine the exon/intron distribution of BnaMYB-CCs through the comparison of CDS and genomic sequences (Guo et al., 2007). The conserved motifs of BnaMYB-CCs were identified using MEME (Bailey et al., 2006) with the flowing parameters: the maximum number of motif: 10; the optimum width: 6–250; the number of repetitions: any. InterProScan was used to annotate the motifs (Quevillon et al., 2005).

Promoter elements analysis

From BRAD, a sequence of 2,000 bp upstream of the start codon was obtained for each MYB-CC gene and submitted to the PlantCARE website to determine the distribution of cis-elements in MYB-CC promoters (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002).

Gene expression patterns in different tissues

The expression patterns of MYB-CC genes in different B. rapa, B. oleracea, and B. napus tissues were analyzed with previously published transcriptomic data (Tong et al., 2013; Yu et al., 2014; Sun et al., 2017). RNA-seq data were collected from the GEO database using accession numbers GSE43245 and GSE42891.

Identification of MYB-CC genes in B. napus, B. rapa, and B. oleracea

The B. napus, B. rapa, and B. oleracea genomes were searched using the protein sequences of 15 MYB-CCs in A. thaliana. A total of 55 MYB-CC genes were identified in B. napus and were equally distributed in the A_n and C_n subgenomes (27 genes in A_n and 28 genes in C_n, respectively) (Table 1). A nomenclature system was used to distinguish the MYB-CC genes. The identified MYB-CC genes in B. napus were designated as BnaMYB-CC01 to BnaMYB-CC55 based on their order on A_n01, C_n01, A_n02, C_n02, etc.. The majority of BnaMYB-CCs contained 148 to 526 amino acids with molecular weights ranging from 16.991 kDa to 59.691 kDa. The predicted isoelectric points (pI) varied from 4.99 to 10.55 (Table 1). A total of 30 MYB-CC genes were identified in B. rapa and were designated as BraMYB-CC01 to BraMYB-CC30 based on their order on A_r01-Ann (Table S1). Furthermore, 34 MYB-CC genes were identified in B. oleracea and were assigned as BolMYB-CC01 to BolMYB-CC34 based on their order on C_o01-Cnn (Table S2).

Evolutionary analysis of MYB-CC family

To explore the evolutionary relationship of the MYB-CC gene family, an unrooted phylogenetic tree was generated using the Neiboring-Joining (NJ) method based on the full-length of 134 MYB-CC protein sequences from A. thaliana, B. napus, B. rapa, and B. oleracea. As shown in Fig. 1, the MYB-CC proteins were classified into nine distinct groups based on the branch of the phylogenetic tree. Group E and Group F each had 20 members, which were the larger groups. Group I was the smallest group with only two MYB-CC members. Group A, B, C, D, G, and H contained 19, 14, 14, 16, 13, and 14 members, respectively. One of the Arabidopsis MYB-CC members (coded by At2g01060) was not found to have a homolog in B. napus. The MYB-CC proteins in Group E were homologous to AtPHR1 and AtPHL4, and the members of Group B were homologs of AtPHL2 and AtPHL3. AtPHL1 belonged to Group F (Fig. 1).

Chromosomal distribution of MYB-CC genes

In order to determine the chromosomal distribution of MYB-CC genes, we searched the Brassica genome database and mapped MYB-CC genes to the corresponding chromosomes. The results showed that 50 BnaMYB-CC genes were distributed across 17 chromosomes except for A_n04 and C_n04 (Fig. 2). There were 27 BraMYB-CC genes located on nine chromosomes except for A_r04, and 28 BolMYB-CC genes were distributed on eight chromosomes except for C_o04 (Fig. 2). Five BnaMYB-CC genes (BnaMYB-CC51, BnaMYB-CC52, BnaMYB-CC53, BnaMYB-CC54, and BnaMYB-CC55), three BraMYB-CC genes (BraMYB-CC28, BraMYB-CC29, and BraMYB-CC30), and six BolMYB-CC genes (BolMYB-CC29, BolMYB-CC30, BolMYB-CC31, BolMYB-CC32, BolMYB-CC33, and BolMYB-CC34) were mapped onto unanchored scaffolds (http://www.genoscope.cns.fr and http://brassicadb.org). In the A_n subgenome of B. napus, chromosome A_n07 carried the most (seven) MYB-CC genes. Chromosome A_n06 and A_n08 only contained one MYB-CC gene and chromosome A_n04 did not have MYB-CC genes. Comparatively, BnaMYB-CC genes were more prevalent in the C_n subgenome. Except for chromosome C_n04, the other chromosomes in the C_n subgenome contained two to four MYB-CC genes (Table S3).

The relative position of seven MYB-CC genes changed on the chromosomes after allotetraploidy. By evolving from diploid to allotetraploid, two BraMYB-CC genes were rearranged in the A_n subgenome belonging to B. napus. Comparatively, the positions of five MYB-CC genes in the C_o and C_n were changed, of which three genes were on chromosome C06 (Fig. 2; Table S4). MYB-CC genes were located on several chromosomes, namely A_r03-A_n03, A_r09-A_n09, A_r10-A_n10, and C_o08-C_n08. This was a conservative process compared to the gene distribution of B. napus and its diploid ancestors (Fig. 2). The results indicated that the genetic variation took place during the evolution of the B. napus genome from its diploid progenitors. The A_n subgenome was more stable than the C_n subgenome, which may be due to the more abundant homologous exchanges or richer transposable elements in the C subgenome (Chalhoub et al., 2014).

Collinearity analysis and gene duplication

Collinearity analysis of the MYB-CC genes was performed on B. napus and its two parental species (B. rapa and B. oleracea). As shown in Fig. 3, the collinear MYB-CC genes were widely distributed in the two subgenomes (A_n and C_n) and two genomes (A_r and C_o), indicating that they promote evolution in the MYB-CC gene family. We also revealed the extension mechanism of the MYB-CC gene family in B. napus by studying gene duplication events. The results showed that 35 BnaMYB-CC genes formed 36 gene pairs with a high homology between gene and protein sequences. Some genes were involved in the duplication events more than once (Table 2). Among these segmental duplication events, 34 events happened between diverse chromosomes, and two events occurred on the same chromosome (BnaMYB-CC26 and BnaMYB-CC29 on C_n06; BnaMYB-CC32 and BnaMYB-CC35 on A_n07) (Fig. 2). Tandem duplication genes were described as a cluster of duplicated genes within 200 kb (Zhu et al., 2020). However, the two pairs of duplicated genes on the same chromosome were not closely related (Fig. 2; Table 1), suggesting that they were not the product of tandem duplication events.

The ka (non-synonymous mutation rate) and ks (synonymous mutation rate) were calculated using the DnaSP program to better understand the effects of evolutionary constraints on BnaMYB-CC genes in B. napus. The ratios of ka and ks of two gene pairs (BnaMYB-CC02 and BnaMYB-CC06; BnaMYB-CC10 and BnaMYB-CC13) were higher than one, indicating that the evolutionary process was actively selected and these genes may function redundantly, resulting in fast evolution (Table 2). However, the remaining 34 gene pairs’ ka and ks rations were less than one, indicating that most BnaMYB-CC genes were conservative and selected in the evolutionary process, resulting in a slower evolutionary rate (Table 2).

Gene structure and conserved motifs

To better understand the diversity of MYB-CC members, the gene exon/intron structure and conserved protein motifs were investigated. All BnaMYB-CC genes were interrupted by several introns, and the exon/intron distribution was complicated (Fig. 4A). BnaMYB-CC37 contained the most introns (11), while several members in group A had the fewest introns (BnaMYB-CC05, BnaMYB-CC09, and BnaMYB-CC22 with four introns each). The members belonging to the same group showed similar exon/intron arrangements (Fig. 4A). For example, the BnaMYB-CC genes in Group H had five introns with similar positions and lengths, suggesting they may have resulted from gene duplication.

Forty-five pairs of homologous genes, which have the closest genetic distance, were analyzed to further investigate the variation of gene structure after allotetraplpidy. The results indicated that the genes were clustered together at the terminal level of the phylogenetic tree (Table S5). Seven pairs of homologous genes showed intron loss/gain variations, while the remaining 38 pairs had the same number of introns (Table S5). In contrast to their ancestral genes, most BnaMYB-CC genes had an identical exon/intron phase (Fig. 4A). It should be noted that intron acquisition events mainly occurred in group C. Specifically, BnaMYB-CC01, BnaMYB-CC20, and BnaMYB-CC54 had one more intron than their homologous genes in the diploid progenitors (Table S5). The majority of the homologous pairs retained consistent number and a similar phase of exon/intron, indicating that MYB-CC gene structures was conserved.

Ten conserved motifs (motifs 1–10) with lengths ranging from 15–50 amino acids were identified using the MEME program (Fig. 4B). The motifs were annotated using the InterProScan program. The results showed that the MYB-like DNA-binding domain and the predicted coiled-coil domain (motif 1 and 2) were found in all BnaMYB-CCs (Fig. 4C). Moreover, similar to the exon/intron organization, the members belonging to the same group also showed similar motif composition. Motif 6 and 7 were found uniquely in the BnaMYB-CC proteins in Group D. Only the proteins of Group A contained motif 8, while motif 10 was found mainly in Group C (Fig. 4C). The other six motifs were widely distributed among the groups. For example, motif 9 was absent in Group A and Group D (Fig. 4C) indicating that there were similarities in the conserved sequences of different phylogenetic groups.

Cis-elements in promoters of MYB-CC genes

To investigate the cis-elements in the promoters of MYB-CC genes, a 2 kb sequence upstream of the start codon of each gene was extracted from the B. napus genome. The promoter sequences were submitted to PlantCARE. As shown in Fig. 5 and Table S6, two types of cis-elements (phytohormone-related and abiotic stress-related elements) were selected. Methyl jasmonate (MeJA)-responsive element (TGACG) and ABA-responsive elements (ABRE) were present in 78.1% and 74.5% of promoters of BnaMYB-CCs, respectively. In the promoters of 55 BnaMYB-CCs, 21 contained auxin-responsive elements (TGA-element, AuxRE, and AuxRR-core), 26 contained SA responsive element (TCA element), and 27 contained gibberellin-responsive elements (GARE-motif, TATC-box, and P-box). Two cis-lements (TC-rich repeats and LTR, related to abiotic stress) were widely distributed in BnaMYB-CC promoters. BnaMYB-CC20 had seven LTRs in its promoter (Fig. 5), indicating that BnaMYB-CC20 may play a vital role in response to low-temperature stresses. The results suggest BnaMYB-CCs may be involved in stress resistance and the hormone signaling pathway. It is worth noting that the cis-elements related to light response were found in significant quantity in the promoters of BnaMYB-CC. However, our work focused on the cis-elements related to phytohormones and thus, light responsive cis-elements are not shown in Fig. 5.

Expression profiles of MYB-CC genes in different tissues

To further understand the biological functions of MYB-CCs, the expression profiles of the genes in different tissues of B. napus, B. rapa, and B. oleracea were assayed based on the transcriptomic data. As shown in Fig. 6, MYB-CC genes were expressed widely in diverse tissues. There was a high accumulation of BraMYB-CCs transcripts found in the roots (Fig. 6A; Table S7). BraMYB-CCs in Group A (BraMYB-CC06, BraMYB-CC17, BraMYB-CC20, and BraMYB-CC28) and Group C (BraMYB-CC03, BraMYB-CC12, and BraMYB-CC26) were extensively expressed in tissues except for silique (Fig. 6A). The transcripts of BolMYB-CCs were generally found in callus (Fig. 6B; Table S7). In B. napus, BnaMYB-CC05, BnaMYB-CC06, BnaMYB-CC19, BnaMYB-CC40, and BnaMYB-CC42 were significantly expressed in all tissues (Fig. 6C). The majority of BnaMYB-CC genes in Group A (BnaMYB-CC09, BnaMYB-CC26, BnaMYB-CC29, BnaMYB-CC32, and BnaMYB-CC35) and Group C (BnaMYB-CC01, BnaMYB-CC20, BnaMYB-CC46, BnaMYB-CC49, and BnaMYB-CC54) showed a low level of expression in the ovule, while their transcripts were found in other tissues (Fig. 6C). BnaMYB-CC12, BnaMYB-CC15, BnaMYB-CC27, BnaMYB-CC33, BnaMYB-CC48, and BnaMYB-CC55 were expressed in fewer tissues (Fig. 6C).

Expression analysis of MYB-CC genes under low Pi stresses

The expression profiles of MYB-CCs under low Pi stress were analyzed by RNA-Seq data to identify potential MYB-CCs with regulatory roles in Pi starvation responses. The original FPKM values of the transcriptome are shown in Table S8. A heatmap was constructed to display diverse expression levels in the roots and shoots under low or high Pi conditions (Fig. 7). The expression of most MYB-CC genes was more active in the roots than that in the shoots, especially BraMYB-CCs and BolMYB-CCs (Fig. 7). The expression patterns of several genes were altered after allopolyploidy. For example, BraMYB-CC24, BolMYB-CC16, and BolMYB-CC31 were mainly expressed in the roots, while their orthologous genes BnaMYB-CC44, BnaMYB-CC26, and BnaMYB-CC23 were expressed dramatically in the shoots (Fig. 7). However, other genes retained the same expression pattern after allopolyploidy, implying that the BnaMYB-CC genes were conserved and maintained similar functions in diploid progenitors (Fig. 7; Table S8).

Approximately half of the BnaMYB-CCs were induced by low Pi stress in both roots and shoots (Fig. 7). Notably, some BnaMYB-CCs were induced twofold or more in roots (BnaMYB-CC01, BnaMYB-CC32, and BnaMYB-CC44) or shoots (BnaMYB-CC04, BnaMYB-CC10, BnaMYB-CC16, and BnaMYB-CC18) under low Pi conditions (Table S8). In B. rapa and B. oleracea, approximately one-third of BraMYB-CCs and two-thirds of BolMYB-CCs genes were up-regulated under low Pi stress in roots (Fig. 7; Table S8). In B. napus, BnaMYB-CCs of Group D (BnaMYB-CC12, BnaMYB-CC15, BnaMYB-CC48, BnaMYB-CC50, BnaMYB-CC51, and BnaMYB-CC53) were significantly down-regulated in both roots and shoots under low Pi stress (Fig. 7C; Table S8). Seven BnaMYB-CCs of Group A were moderately induced by low Pi stress in roots, while their expression levels increased in shoots under high and low Pi conditions (Fig. 7; Table S8). However, BnaMYB-CC10, BnaMYB-CC13, BnaMYB-CC27, and BnaMYB-CC34 of Group H were induced by low Pi stress in roots prominently. All six members of BnaMYB-CCs in Group H were induced by low Pi stress in shoots (Fig. 7; Table S8). The results indicate that the expression patterns were conservative for MYB-CCs of the same evolutionary group.