The B. rapa annotated genome sequences were downloaded from the Brassica database (BRAD; http://brassicadb.org/brad/; Cheng et al., 2011). Several datasets and multiple steps were used to conduct a whole-genome search. Hmmsearch from the HMMER program (3.0) (Eddy, 1998) was used to search the Brassica rapa genome to identify putative ABCs containing the ABC-transporter domain (PF00005), the ABC-2 transporter domain (PF01061), the ABC transporter transmembrane region domain (PF00664), the CYT domain (PF01458) or the mce related protein domain (PF02470) at a threshold of E < 1E⁻⁴. The Hidden Markov Model (HMM) profiles of these domains were downloaded from the Pfam 27.0 database (http://Pfam.sanger.ac.uk/; Finn et al., 2013). Then, we manually verified the candidate ABCs using the SMART database (http://smart.embl-heidelberg.de/; Letunic et al., 2012), and sequences with obvious errors and/or lengths of less than 100 amino acids were manually removed.

Amino acid sequences of the ABC-transporter domain (PF00005) were extracted according to their locations in ABC proteins. Multiple alignments of ATPase domain sequences were conducted using ClustalW (Thompson et al., 1997) of the MEGA6 (Tamura et al., 2013) software with the default options. The phylogenetic trees were constructed using the maximum likelihood (ML) method in MEGA6 with the Jones, Taylor, and Thornton (JTT) amino acid substitution model. One thousand bootstrap replications were performed in each analysis. The remaining parameters were set to the defaults.

Conserved motifs of ABC proteins were identified using the local MEME software version 4.9.0 (Bailey et al., 2009), with the following parameter settings: maximum number of motifs = 15, and optimum motif width = 10–200 residues; the other parameters were the defaults. The coding sequences (CDSs) and the DNA sequences were used to predict the exon-intron structure of ABCs using the Gene Structure Display Server 2.0 (GSDS 2.0) (http://gsds.cbi.pku.edu.cn; Hu et al., 2015).

All of the promoter sequences (2,000 bp upstream of initiation codon “ATG”) of BraABCs were extracted from the B. rapa genome according to the Generic File Format (GFF) file. Then, the cis-regulatory elements of promoters for each gene were identified by PLACE Web Signal Scan-PLACE (https://sogo.dna.affrc.go.jp/cgi-bin/sogo.cgi?lang=en&pj=640&action=page&page=newplace). The interaction network of BraABC proteins was constructed by applying STRING software (Search Tool for the Retrieval of Interacting Genes/Proteins, http://string-db.org/).

Chromosomal locations of BraABCs were drawn from top to bottom on B. rapa chromosomes according to their positions in genome annotation. The Perl in-house program was used to draw the location images of BraABCs on chromosomes. Syntenic relationships of ABCs between Arabidopsis and B. rapa were identified by searching “syntenic gene” in the Brassica database (BRAD) (Cheng et al., 2011). Then, the Circos program (Krzywinski et al., 2009) was applied to present the syntenic relationships on their chromosomes.

After the identification of the BraABCs, syntenic BraABCs paralogs pairs among the three B. rapa subgenomes were identified by searching “syntenic gene” in BRAD (Cheng et al., 2011). The non-synonymous substitution rate (Ka), the synonymous substitution rate (Ks) and the ω (= Ka/Ks) of paralog pairs were estimated by KaKs_Calculator 2.0 (Wang et al., 2010). The duplication date of paralog pairs was estimated according to formula T = Ks/2λ, assuming a clock-like rate (λ) of 1.5 synonymous substitutions per 10⁸ years for B. rapa (Koch et al., 2000).

To analyze the expression patterns of BraABCs in different tissues, we used the reported Illumina RNA-seq data by Tong et al. (2013). Of the dataset, five tissues, including the roots, stems, leaves, flowers and siliques of B. rapa accession Chiifu-401-42 were analyzed. Gene expression levels were quantified by FPKM (Fragments Per Kilobase of transcript per Million fragments mapped) values. An expression cluster of BraABCs in the five tissues was generated by Cluster 3.0 (http://bonsai.hgc.jp/~mdehoon/software/cluster/software.htm), and the cluster results were shown using Tree View software (http://jtreeview.sourceforge.net/). The RNA-seq dataset was also used to analyze the expression profiles of syntenic paralog pairs in different tissues.

To investigate the expression patterns of syntenic BraABC paralog pairs under normal and abiotic stress conditions, the cultivar Chiifu-401-42 was cultivated in potting soil at 24/16°C, with a photoperiod of 14/10 h for day/night in a growth chamber. Five-leaf-stage plants were used to conduct different abiotic stress treatments under continuous time (6, 12, and 24 h). Hormone treatment was performed with 100 μM ABA; for drought stress treatment, the pots were irrigated with 15% (w/v) polyethylene glycol PEG and kept standing in the irrigation solution for 30 min under normal growth conditions; salt treatment was performed with 250 mM NaCl. After the treatments, the total RNAs were isolated from young leaves using an RNA kit (TaKaRa, Dalian, China). The RNAs were reverse transcribed into cDNA with a PrimeScript RT reagent Kit (TaKaRa, Dalian, China). Six pairs of syntenic paralogs were randomly selected to conduct qRT-PCR. The B. rapa actin gene Bra028615 was used as an internal control. Primers for qRT-PCR were designed using Beacon Designer 7. The qRT-PCR experiments were performed with three biological and technical replicates. The SYBR Premix Ex Taq kit® (TaKaRa, Dalian, China) was used to detect gene expression according to the manufacturer's recommendations on the One-step Real-Time PCR System Time PCR Detection System (Applied Biosystems, Foster City, CA, USA). The cycling profile was as follows: 94°C for 30 s, 94°C for 10 s for 40 cycles, and 58°C for 30 s, followed by a melting curve analysis at 65°C for 10 s for 61 cycles. The relative gene expression levels were analyzed using the comparative Ct value method (Heid et al., 1996).

To explore all of the putative ABCs in the B. rapa genome, the HMM profiles of the ABC-transporter domain, the ABC-2 transporters domain, the ABC transporter transmembrane region domain, the CYT domain and the mce related protein domain were employed as queries to search against the protein database of B. rapa using the Hmmsearch program. A total of 179 BraABCs, which contained at least one of the above domains, were identified (Table 1). Among the 179 BraABC proteins, 162 proteins contained one or more NBD domains, and an additional 17 proteins without NBD domains were assigned to ABCI, based on the domain composition of the ABC proteins (Verrier et al., 2008).

Understanding the evolutionary relationship between BraABC proteins and AtABC proteins could contribute to classifying and assessing the potential functions of BraABCs. The NBD domains of the 162 BraABC proteins and 117 AtABC proteins (117 AtABC proteins contained the NBD domain) were aligned, and an unrooted phylogenetic tree was constructed by MEGA6 using the Maximum Likelihood (ML) method with 1,000 bootstraps. The AtABC proteins of each subfamily were grouped together, whereas the proteins of ABCI were dispersive in the phylogenetic tree. Overall, BraABC proteins followed the phylogenetic pattern of Arabidopsis (Supplementary Figure 1). Based on the phylogenetic relationship with AtABCs, BraABCs were also divided into eight subfamilies: ABCA, ABCB, ABCC, ABCD, ABCE, ABCF, ABCG and ABCI. ABCG, ABCB, ABCI, and ABCC were the largest subfamilies containing 63, 38, 30, and 21 members, respectively. ABCA, ABCE, and ABCF included 11, 7, and 7 members, respectively. Only 2 members were identified in ABCD (Table 1). The protein domain organizations of BraABCs in each subfamily were consistent with that of Arabidopsis (Verrier et al., 2008).

Similarly, Hmmsearch searches were also conducted in the selected angiosperms: P. trichocarpa, G. max, A. lyrate, C. papaya, V. vinifera, B. distachyon, O. sativa, and A. trichopoda. The number of identified ABCs in these species ranged from 113 to 271 (Figure 1). There were 271 ABCs existing in the G. max genome, while the members in other species ranged from 113 in C. papaya, 130 in Arabidopsis, 132 in A. lyrate, 137 in A. trichopoda, 138 in B. distachyon, 141 in O. sativa, 179 in B. rapa, 181 in V. vinifera and 204 in P. trichocarpa. The identified ABCs in these species are listed in Supplementary Table 1.

To gain insight into the evolutionary history and relationship of ABCs in plants, we constructed a phylogenetic tree using ABCs of six representative species, including Arabidopsis (Verrier et al., 2008), B. rapa, C. papaya, V. vinifera, O. sativa, and A. trichopoda. In contrast with the ABCA-ABCG members, some ABCI proteins lacked the NBD domain, and some ABCI members were scattered in the phylogenetic tree, although they contained the NBD domain. Therefore, our phylogenetic tree was constructed using the NBD domains of ABCA-ABCG members from these species. As shown in Figure 2, seven distinct clades were clearly distinguished based on the tree topology, namely ABCA, ABCB, ABCC, ABCD, ABCE, ABCF, and ABCG, suggesting that these groups have diverged prior to the divergence of angiosperm. Next, the copy number of ABCs in each subfamily of these species was quantified (Figure 1). In each studied species, the copy numbers of different subfamilies varied greatly. For the copy numbers of subfamilies in all of the studied species, we found that ABCB, ABCC, ABCG, and ABCI were the larger subfamilies. ABCA, ABCD, ABCE, and ABCF contained a small number of members, and the copy numbers of ABCD, ABCE, and ABCF were relatively constant. B. rapa has undergone a further WGT since its divergence from Arabidopsis. Compared with Arabidopsis, the expansion of the ABC gene family in B. rapa was largely the result of the expansion of ABCs in ABCB, ABCC, ABCG, and ABCI (Figure 1).

The 179 BraABCs were designated as BraABCA1 through BraABCI30 according to their respective subfamily and their location on the chromosome (Table 1). The comprehensive information of BraABCs, including the locus name, gene location, predicted protein length, exon number, molecular weight (MW), isoelectric point (PI) and transmembrane helices number (TMs), is listed in Table 1. The size of the BraABC proteins varied from 118 (BraABCB31) to 5,408 (BraABCI6) amino acids, with an average of ~984 amino acids. There were few differences in the protein lengths of ABCC, ABCE and ABCF, whereas the protein lengths of ABCB, ABCG and ABCI changed greatly. The MWs of the BraABC proteins ranged from 12.87 kDa (BraABCB31) to 611.88 kDa (BraABCI6). The variation range of PI was from 4.79 (BraABCI18) to 10.52 (BraABCI22), indicating that these BraABC proteins might exist and function in different parts of the cells. Most of the BraABC proteins (162 out of 179, 90.5%) presented multiple exons and sixty proteins (33.5%) that contained more than 10 exons. The majority of the BraABC proteins contained transmembrane helices (TMs), except for all of the ABCE and ABCF proteins as well as parts of the ABCI proteins. The ubiquity of TMs was in accordance with the transmembrane transport functions of ABCs. Pairwise comparisons of the full-length BraABC protein sequences within each subfamily were conducted. The protein sequences within ABCA, ABCB, and ABCI exhibited higher identities than other subfamilies (Supplementary Figure 2). In contrast, the sequence identities of ABCD and ABCG were lower, indicating that the degree of sequence divergence of the two subfamilies was higher than other subfamilies. The other families, such as ABCC, ABCE, and ABCF exhibited medium sequence identities.

A total of 173 BraABCs were located on the ten chromosomes (96.7%), and the other 6 genes were located on unanchored scaffolds (Figure 3). The BraABCs were unevenly distributed on the ten chromosomes. The largest number of BraABCs was found on chromosome A09 (26 genes), followed by A03 (25 genes). Chromosomes A06 and A05 also contained more than 20 BraABCs (22 and 21, respectively). In contrast, chromosomes A04, A01, A07, A10, A08, and A02 contained 18, 14, 14, 12, 11, and 10 BraABCs, respectively. Our results indicated that BraABCs within each subfamily were also unevenly distributed on the ten chromosomes.

Simultaneously, we further identified tandem arrays of BraABCs along the 10 B. rapa chromosomes. A tandem array was defined as multiple members of BraABCs occurring within the same or neighboring intergenic regions. As shown in Figure 3, 11 BraABCs clusters (genes marked in purple) containing 28 tandemly duplicated genes were identified on chromosomes A03, A04, A05, A06, and A09.

To better understand the conservation and diversity of motif compositions and gene structures of BraABCs, the conserved motifs and exon-intron organization of BraABCs were analyzed.

The conserved motifs of BraABC proteins in ABCA-ABCG were analyzed using MEME software, and we identified 15 conserved motifs (Figure 4). The lengths of the conserved motifs were between 26 and 154 amino acids. The number of the conserved motifs in each BraABC protein varied from 1 to 8. Several motifs were widely spread among BraABC proteins; for instance, motifs 1 and 3 were present in the most proteins of ABCA-ABCG. In contrast, other motifs were specific to only one or two subfamilies. For example, motifs 4, 5, 8, and 12 were specific to ABCB and ABCC, motifs 10 and 15 were specific to ABCB, and motif 14 was specific to ABCA. ABCG contained the specific motifs 2, 6, 7, 9, 11, and 13, and these motifs were probably required for specific protein functions. The varieties of motif compositions in different subfamilies suggested sources of functional differentiation in BraABCs during the evolutionary process. BraABC proteins grouped into the same subfamily exhibited similar motif characteristics, suggesting functional similarities for members in the same subfamily. Logos of the 15 identified motifs are shown in Supplementary Figure 3.

We also analyzed the exon-intron organizations of BraABCs using GSDS 2.0. Similar to the motif composition, most of the BraABCs within the same subfamily generally exhibited similar gene structure in terms of the numbers and lengths of introns and exons. Moreover, BraABCs that had a closely phylogenetic relationship and shared a more similar gene structure (Supplementary Figure 4).

To explore the influences of WGT on the evolution of BraABCs, we studied the conservation of BraABCs after WGT. After diverging from Arabidopsis, the gene number of the B. rapa genome was ~42,000, which is considerably lower than the simple triplication of the ~30,000 genes in the Arabidopsis genome. Therefore, substantial genes must have been lost (fractionation) in B. rapa after triploidization (Wang et al., 2011). To evaluate the fractionation extent of BraABCs, we compared the retention of BraABCs relative to the set of 3,580 neighboring genes (10 on either side), flanking the 179 BraABCs (Supplementary Table 2). In comparison, 47.2% of BraABCs were retained in one copy, which was greater than the retention of neighboring genes (38.7%). A similar percent (33.7 and 30.6%, respectively) of BraABCs and neighboring genes were retained in the two copies. And more BraABCs (19.1%) were retained in the three copies than their neighboring genes (9.2%) (Figure 5A).

The triplicated B. rapa genome can be divided into three subgenomes: LF, MF1, and MF2. The subgenomes were grouped based on the level of gene loss relative to Arabidopsis since the WGT 13–17 MYA. Compared with their neighboring genes, more BraABCs were significantly retained in all three subgenomes. Consistent with the whole genome level, BraABCs distributed in the LF subgenome were more than two other subgenomes (Figure 5B). Given the well-conserved synteny between Arabidopsis and B. rapa (Cheng et al., 2012), in our study, there were four collinear relations for ABCs between B. rapa and Arabidopsis: a single copy of BraABC to a single copy of Arabidopsis (SB to SA); a tandem array of BraABCs to a tandem array of Arabidopsis (TB to TA); a single copy of BraABC to a tandem array of Arabidopsis (SB to TA); and a tandem array of BraABCs to a single copy of Arabidopsis (TB to SA). The four syntenic relations were presented using the Circos program (Figure 5C). On the whole, BraABCs were preferentially retained, consistent with the gene balance hypothesis and demonstrating the significance of ABCs in B. rapa.

To identify cis-elements related to transcriptional regulation, 2 kb promoter regions for each of the 179 BraABCs were identified. Gaps existed in the promoter regions of 13 BraABCs after assembling the B. rapa genome. Therefore, these BraABCs were excluded in the following study. Next, we applied the Plant Cis-acting Regulatory DNA Elements (PLACE) website to analyze the cis-elements and identified a total of 291 different cis-elements in all of the studied BraABCs.

A total of 13 common cis-regulatory elements were identified in all of the promoter regions of the BraABCs, which were highly conserved among all of the studied BraABCs (Table 2). Two common cis-regulatory elements, WRKY71OS and GT1CONSENSUS were responsive to plant hormones, including ABA, GA, and SA, suggesting that ABA, GA, and SA could affect the expression levels of BraABCs. Some common cis-regulatory elements were responsive to stresses, such as pathogens (WRKY71OS) and wounds (WBOXNTERF3). Of the 13 common cis-regulatory elements, GATABOX, IBOXCORE, and GT1CONSENSUS were required for transcriptional regulation by light. DOFCOREZM participated in carbon metabolism, and GATABOX played a role in molecular light switching; therefore, we speculated that BraABCs were likely to participate in energy metabolism. Two common cis-elements, POLLEN1LELAT52 and GTGANTG10, were required for transcriptional regulation in pollen, suggesting that BraABCs might be involved in the reproductive process.

We also identified unique cis-element sequences that are present in the promoter regions of unique BraABCs. In total, 40 unique cis-elements were identified in the promoter regions of the 36 BraABCs (Supplementary Table 3). Interestingly, 3 unique cis-elements (ACGTSEED3, ABREMOTIFIOSRAB16B, and VOZATVPP) were presented in only BraABCF3. BraABCG18 and BraABCB10 both contained two unique cis-elements. There was only one unique cis-element present in the other 33 BraABCs. The unique cis-element of these genes indicated their expression specificity.

To characterize the expression patterns of individual BraABCs in different tissues, we used publicly available RNA-Seq data of different tissues as a resource (Tong et al., 2013). A heat map displaying the expression profiles of BraABCs in roots, stems, leaves, flowers and siliques was generated, and the BraABCs were clustered by their expression patterns. As shown in Figure 6A, the majority of BraABCs presented different expression patterns, whereas a few of them exhibited similar expression patterns. A total of 165 BraABCs (92.2%) were determined as being expressed in at least one tissue, and 117 BraABCs (65.4%) expressed in all tissues, including 3 ABCA genes, 26 ABCB genes, 19 ABCC genes, 2 ABCD genes, 5 ABCE genes, 6 ABCF genes, 37 ABCG genes and 19 ABCI genes. Relatively, the expression of 14 BraABCs was not detected in any of the five tissues. Some BraABCs also exhibited tissue-specific expression, and there were 8, 4, 4, 5, and 1 tissue-specific BraABCs in the roots, stems, leaves, flowers and siliques, respectively (Figure 6B), suggesting that these genes may play specific roles in the relevant tissues. As shown in Figure 6A, some genes of ABCA, ABCB, and ABCG exhibited variable expression profiles throughout the five tissues, such as BraABCA10, BraABCB13, and BraABCG25. In contrast, the majority of genes in ABCC, ABCD, ABCE, ABCF, and ABCI presented relatively consistent expression levels across the five tissues, such as BraABCC18, BraABCD2, BraABCE6, BraABCF1, and BraABCI14. As shown in Figure 6A, we found the expression profiles of BraABCs in the roots exhibited large differences compared with the other four tissues, whereas similar expression patterns of BraABCs were observed for the stems and leaves, and the flowers and siliques.

The three subgenomes of B. rapa are the evolutionary products of WGT and a lot of synteny blocks among them. Syntenic paralogs are genes that are located in the syntenic fragments. Syntenic BraABC paralog pairs among LF, MF1, and MF2 were identified by searching “syntenic gene” in the BRAD (Cheng et al., 2011). In total, 76 pairs of BraABC syntenic paralogs were identified among the three subgenomes (Table 3). Synonymous (Ks) and non-synonymous (Ka) values were calculated to explore the selective pressures on these paralog pairs. We found that the ω (= Ka/Ks) ratios of 10 syntenic paralogs (13.16%) were more than 1, indicating positive selection on these BraABC paralogs. The ABCC21_ABCC3 pair had the highest ω ratios, with a ratio of 1.9264. The rest of the paralogs were under purifying selection, with ω ratios lower than 1. Among them, the ω ratios of 22 paralog pairs were lower than 0.1, suggesting that these gene pairs underwent strong purifying selection stress, making the functions of these paralog pairs trend toward relative similarity. Furthermore, the duplication time of these paralogs pairs were estimated by using a relative Ks measure as a proxy for time. The estimated duplication time of BraABC paralog pairs indicated that it spanned from 3.45 to 42.2 MYA, with an average duplication time of ~15.13 MYA, which was consistent with the recent WGT time of B. rapa (13–17 MYA) (Wang et al., 2011).

Profiling gene expression may provide an ample amount of information about the mode and tempo of duplicated genes. To predict the evolutionary fates of syntenic BraABC paralog pairs, we took advantage of the deep RNA-Seq data of B. rapa to investigate their expression profiles in different tissues, including roots, stems, leaves, flowers and siliques (Tong et al., 2013). Generally, if a Pearson correlation coefficient (PCC) of the two factors is greater than 0.6, we can infer that they have a close positive correlation. Among the 76 pairs of paralogs, 27 pairs had positive expression correlations across the five tissues, which were greater than 0.6 measured by the PCCs (Figure 7). Such high positive correlations among these paralog pairs probably indicated functional conservation or sub-functionalization for these paralog pairs after duplication. Simultaneously, 40 pairs of syntenic paralogs were negatively correlated with negative PCC values or exhibited little correlation with small positive PCC values in gene expression, suggesting neo-functionalization for these paralog pairs. Furthermore, there were 9 pairs of paralogs in which one copy of each pair did not express in any of the 5 tissues. Therefore, there were no available (NA) results of PCCs for these paralog pairs. Among the 9 paralog pairs, 8 pairs were composed of genes from two different subfamilies, and one copy of each of these paralog pairs was the ABCI gene. Interestingly, all of these ABCI genes were not expressed in any of the studied tissues, indicating that they probably became pseudogenes. Therefore, pseudogenization of these 9 pairs of paralogs may occur during evolution.

From the cluster of different tissues, the leaves clustered with the stems and the flowers clustered with the siliques in the heat map of BraABC paralog pairs indicated a similar expression profile (Figure 7). Taking an overall view of the cluster for BraABC paralog pairs, paralog pairs with similar evolutionary fates tended to further cluster together, confirming that the expression patterns of paralog pairs were related to their evolutionary fates after duplication. Approximately 42.1% (32 out of 76) of the paralog pairs presented no less than a two-fold expression difference (absolute value of log₂ FPKM ratio ≥ 1) in all of the studied tissues. In brief, the expression correlation analysis of syntenic paralog pairs showed their functional roles in functional conservation, sub-functionalization, neo-functionalization and pseudogenization in the BraABC gene family.

We also investigated the expression patterns of syntenic BraABC paralog pairs under ABA, drought and salt stresses with qRT-PCR. Three pairs of conserved or sub-functional paralogs and three pairs of neo-functional paralogs were randomly selected to study their expression patterns. Primers for qRT-PCR are listed in Supplementary Table 4. As shown in Figure 8, a range of expression profiles of the selected BraABCs were observed after exposure to ABA, drought and salt stress. Some paralog pairs exhibited similar expression patterns, whereas others presented different expression profiles. The paralog pairs with positive expression correlations measured by PCCs across the five tissues also had positive correlations under the stress treatments (Figures 8A–C), which also indicated the functional conservation or sub-functionalization in these paralog pairs. Similarly, paralog pairs which exhibited negative correlations across the studied tissues also presented negative correlations under the stress treatments (Figures 8D–F), suggesting that these paralog pairs were neo-functional. On the whole, the expression correlations of the selected syntenic BraABC paralog pairs across the five tissues were consistent with those under the stress treatments.

The interaction network of BraABC proteins, including the functional and physical interactions, were examined using STRING software and the corresponding database to retrieve the protein interactions. As shown in Figure 9, a dense interaction network formed among the BraABC proteins. Most of the BraABC proteins were involved in interaction relations with other proteins, except for BraABCA4, BraABCA9, BraABCG18, and several ABCI proteins. Interaction relations not only existed among proteins within the same subfamily but also among proteins from different subfamilies. Though interacting frequently with proteins in other subfamilies, some ABCB proteins did not interact with proteins of ABCB, shown as being surrounded by the yellow curve in Figure 9. The majority of proteins located in the center of the network were ABCG proteins, indicating that these proteins had more complex interaction relations with other BraABC proteins. The gene dosage hypothesis predicts that genes will be preferentially retained if their products are dose sensitive, interacting either with other proteins or in networks (Thomas et al., 2006; Birchler and Veitia, 2007). Consequently, ABCG genes were probably more preferentially retained compared with genes of other subfamilies during evolution. Indeed, ABCG had the largest number of genes among all of the BraABC subfamilies.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.