Comparative expression analysis of heavy metal ATPase subfamily genes between Cd-tolerant and Cd-sensitive turnip landraces.

The heavy metal ATPase (HMA) subfamily is mainly involved in heavy metal (HM) tolerance and transport in plants, but an understanding of the definite roles and mechanisms of most HMA members are still limited. In the present study, we identified 14 candidate HMA genes named BrrHMA1-BrrHMA8 from the turnip genome and analyzed the phylogeny, gene structure, chromosome distribution, and conserved domains and motifs of HMAs in turnip (Brassica rapa var. rapa). According to our phylogenetic tree, the BrrHMAs are divided into a Zn/Cd/Co/Pb subclass and Cu/Ag subclass. The BrrHMA members show similar structural characteristics within subclasses. To explore the roles of BrrHMAs in turnip, we compared the gene sequences and expression patterns of the BrrHMA genes between a Cd-tolerant landrace and a Cd-sensitive landrace. Most BrrHMA genes showed similar spatial expression patterns in both Cd-tolerant and Cd-sensitive turnip landraces; some BrrHMA genes, however, were differentially expressed in specific tissue in Cd-tolerant and Cd-sensitive turnip. Specifically, BrrHMA genes in the Zn/Cd/Co/Pb subclass shared the same coding sequence but were differentially expressed in Cd-tolerant and Cd-sensitive turnip landraces under Cd stress. Our findings suggest that the stable expression and up-regulated expression of BrrHMA Zn/Cd/Co/Pb subclass genes under Cd stress may contribute to the higher Cd tolerance of turnip landraces.

Introduction

Heavy metal (HM) contamination is a global environmental problem. Cd, which is highly mobile and extremely toxic, has become the most prominent heavy metal pollutant in soil (Moreno-Caselles et al., 2000, Wei et al., 2006). Plants absorb HMs from soil and transfer them to animals and humans via the food chain (Tauqeer et al., 2016), and are therefore a critical link between HMs and human health. Accordingly, great interest has been focused on decreasing the bioavailability of HMs such as Cd in soils (Buendia-Gonzalez et al., 2010, Sarwar et al., 2017). At the same time, plants capable of HM hyperaccumulation provide opportunities for an efficient, environmentally-friendly method of soil remediation, i.e., phytoremediation (Ehsan et al., 2014, Kramer, 2005). Thus, understanding the genetic mechanisms underlying HM accumulation has become a major goal of plant biologists (Gustin et al., 2009, Zhang et al., 2016). To date, several gene families have been reported to participate in plant absorbing and transporting HMs, such as P-type ATPase, cation diffusion facilitator gene family, Zrt/Irt-like protein, and ATP-binding cassette transporter (Zhang et al., 2013).

The heavy metal ATPases (HMAs), also known as the P1B-type ATPases, play an important role in transporting metals in plants (Li et al., 2015, Zhang et al., 2018b). HMA genes can selectively absorb and transport either metal ions essential for plant growth and development (e.g. Cu2+, Zn2+, and Co2+) or HMs not essential to plants (e.g., Cd2+ and Pb2+) (Li et al., 2015). Typical HMAs contain approximately 6–8 transmembrane helices, a soluble nucleotide binding domain, phosphorylation domain, and a soluble actuator domain (Li et al., 2015). In addition, both sides of the N-terminal and C-terminal metal binding sites include metal binding domains, and P1B-type ATPase is also located in this site, which is a heavy metal-associated regulatory domain as well (Williams and Mills, 2005, Arguello et al., 2007). HMA proteins have been studied at the genomic scale in Arabidopsis thaliana, Oryza sativa, and Populus trichocarpa, which contain 8, 9, and 12 full-length HMA genes, respectively (Williams and Mills, 2005, Li et al., 2015). Based on metal substrate specificity, these HMAs can be classified into two major phylogenetic subclasses: the Zn/Cd/Co/Pb subclass and the Cu/Ag subclass (Axelsen and Palmgren, 1998). AtHMA1–4 in A. thaliana, OsHMA1–3 in O. sativa, and PtHMA1–4 in P. trichocarpa belong to the Zn/Co/Cd/Pb subclass, whereas AtHMA5–8, OsHMA4–9, and PtHMA5.1–8 belong to the Cu/Ag subclass (Williams and Mills, 2005, Li et al., 2015). The functions of HMA members have been analyzed in an increasingly diverse set of species (Wong et al., 2009, Deng et al., 2013, Huang et al., 2016, Zhang et al., 2016). These studies show that HMAs in maintaining homeostasis of trace elements under normal conditions and improving tolerance and accumulation of HMs under HM stress (Wong et al., 2009, Deng et al., 2013, Huang et al., 2016, Zhang et al., 2016).

Brassicaceae is an excellent model for studying the roles of metal ion transporters. Brassicaceae contains species that differ in their tolerance to abiotic stresses such as salinity and trace metals (Assuncao et al., 2003, Prasad and Freitas, 2003, Inan et al., 2004, Megdiche et al., 2007), and a large number of HM hyperaccumulators (Assuncao et al., 2003). Previous studies have extensively examined HM tolerance in different Brassicaceae species (including several varieties of B. campestris, Brassica rapa, Brassica napus, Brassica oleracea and B. carinata), to Zn, Pb, and Cd (Hernandez-Allica et al., 2008). However, among these HMs, Cd accumulation in turnip has attracted more attention (Wei et al., 2006). Turnip (B. rapa var. rapa), a Brassicaceae biennial plant, is widely cultivated in America, Europe, and Asia. To date, few studies have investigated the accumulation characteristics of turnip to various HMs, such as Zn, Cu, Mn, Pb, and Cd (Siddiqui et al., 2014, Parveen et al., 2015, Li et al., 2016). Turnip has been designated a high-Cd-accumulation species (Arthur et al., 2000). Recently, we also found that Chinese turnip has high capacity in absorption and transportation of Cd and some turnip landraces have met the standard for Cd hyperaccumulators (Li et al., 2016, Li et al., 2017a). Moreover, turnip landraces show excellent Cd absorption and extraction abilities compared with three other high-Cd-accumulating species, i.e., B. napus, Phytolacca americana, and Bidens pilosa (Li et al., 2017a). Thus, the genes from turnip involved in Cd accumulation may be important for molecular-assisted breeding of plant resources for phytoremediation. However, the genetic background responsible for differences in Cd accumulation among landraces has not been identified. In the present study, we aimed to identify genes that play a role in Cd hyperaccumulation in B. rapa var. rapa. For this purpose, we systematically analyzed the sequence and structural characteristics of putative HMAs in turnip (Cheng et al. 2016) and we compared the sequence characterization and expression differences of these HMA genes between Cd-tolerant and Cd-sensitive turnip landraces. This study can improve the understanding of the roles of turnip HMAs in Cd tolerance and provide candidate genes for further functional verification or utilization.

Materials and methods

Identification and characteristics analysis of the HMA gene subfamily in turnip

The eight Arabidopsis HMA gene sequences were downloaded from TAIR (http://www.arabidopsis.org) to search turnip HMAs as previously reported method (Yin et al., 2017). All protein sequences containing any of the representative domains of HMA proteins were extracted as candidates. The candidates were then used to search against the GenBank nonredundant protein database. ClustalW was used for sequence alignment between turnip and Arabidopsis, and phylogenetic trees were constructed using the MEGA 7.0 software with the neighbor-joining method (1000 bootstrap test replicates) (Kumar et al., 2016).

The intron/exon distribution diagram of turnip HMA genes was analyzed by the online Gene Structure Display Server (http://gsds.cbi.pku.edu.cn/) (Hu et al., 2015). The chromosomal locations of turnip HMA genes were mapped according to the gene position information by the TBtools software (https://github.com/CJ-Chen/TBtools).

The Pfam tool (http://pfam.xfam.org/search#tabview=tab1) and MEME program (http://meme-suite.org/tools/meme) were used to search for conserved domains and motifs in the turnip HMA protein sequences, respectively (Sonnhammer et al., 1998, Bailey et al., 2006); subsequently, the domain and motif diagrams were drawn using the TBtools software. The putative transmembrane helixes were predicted using TMHMM Server V. 2.0 (http://www.cbs.dtu.dk/services/TMHMM/). The physicochemical parameters of the proteins, including the molecular weight (MW), theoretical isoelectric point (pI), and grand average of hydropathicity (GRAVY) were calculated using the ProtParam tool of ExPaSy (http://web.expasy.org/protparam/) (Artimo et al., 2012).

Cd tolerance identification among turnip landraces

Sixteen Chinese turnip landraces from different regions were selected (Table S1 in Supporting Information) and mature seeds of each landrace were sterilized in 1% NaClO solution for 5 min before washing three times with deionized water. Water on the surface of the seeds was blotted using filter paper. Culture dishes with two layers of filter paper were used as germinating beds. To test Cd tolerance, we used 50 mg L−1 Cd2+ concentration (supplied by CdCl2•2.5H2O); control groups were treated with deionized water (CK). Each culture dish was sown with 20 seeds that received either 5 ml of Cd2+ solution or the water control. There were three duplicates for each group. Cd treatments were performed in an incubation room (22–24 °C, 12-h light/12-h darkness). Culture dishes were sealed to minimize water loss, and water was supplemented as needed to keep the filter paper moist. Seeds with 2 mm radicles were considered to have germinated. The number of germinated seeds was recorded 7 d after germination, at which point ten seedlings in each dish were washed thrice using deionized water and transplanted into Cd-free soils. Two weeks later, the number of surviving seedlings for each sample was recorded. The germination percentage index (GPI) and seedling survival index (SSI) were calculated according to the following formulas:

Germination percentage index = germination percentage of treatment sample/germination percentage of control sample

Seedling survival index = seedling survival percentage of treatment sample/seedling survival percentage of control sample

Plant preparation and treatment

Seeds of selected Cd-tolerant and Cd-sensitive turnip landraces were sown in soil pots (three seeds per pot) under greenhouse condition (22–25 °C, 12-h light/12-h darkness, 40%∼50% relative humidity). After one month the plants were transferred into hydroponic containers. Each of two containers was equipped with 10 L solutions and was planted with 12 Cd-tolerant and 12 Cd-sensitive turnip landraces (24-h oxygenation). The plants were cultured using water for 48 h and subsequently cultured using 1/2 Hoagland's culture medium (Xue et al., 2014) for 72 h. Then, fresh 1/2 Hoagland's culture medium was changed and one group of samples were added with 2.5 mg L−1 Cd2+ (supplied by CdCl2 • 2.5H2O) for 6 h. The roots and leaves of each group were harvested, grinded into powder, and then stored at −80 °C for RNA isolation.

RNA extraction and cDNA synthesis

Total RNA samples were isolated using the Eastep® Super Total RNA Extraction Kit (Promega, Madison, WI, USA). RNA concentration was determined by NanoDrop1000 (NanoDrop Technologies, Inc.) and the integrity was checked on 0.8% agarose gel. A total of 3 μg RNA was reverse transcribed using the GoScript Reverse Transcription System (Promega, Madison, WI, USA) to generate cDNA.

Gene cloning

The forward and reverse primers of Zn-subclass BrrHMA genes were generated according to the CDS sequences (Supplementary Table S2). cDNA of both Cd-tolerant and Cd-sensitive turnip landraces was used as template for each gene. PCR parameters were 95 °C = for 5 min followed by 30 cycles of 95 °C for 30 s, 51–58 °C for 30 s and 72 °C for 2 min, which was then followed by a 10 min extension at 72 °C (Table S2 in Supporting Information). Purified PCR products were ligated to pGEM-T Easy vector and subsequently transformed to E. coli DH5α. Positive clones were sequenced.

Gene expression analysis

Optimal forward and reverse primers of Zn-subclass BrrHMA genes were designed (Table S3 in Supporting Information) through the online tool Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) for qRT-PCR analysis. qRT-PCR was conducted in triplicate with different cDNAs from different tissues and treatments using FastStart Universal SYBR Green Master (Rox, Roche, Indianapolis, IN, USA) and a 7500 Sequence Detection System (Applied Biosystems, USA). The reaction parameters for thermal cycling were as follows: 95 °C for 10 min, followed by 40 cycles of 94 °C for 5 s, and 60 °C for 15 s. Turnip beta-tubulin gene was amplified as an internal control. Relative gene expression levels were obtained by dividing the extrapolated transcript levels of the target genes by the levels of the internal control from the same sample.

Statistical analysis

Independent sample t-test was performed to analyze significant differences between two samples using SPSS version 18.0.

Results

Identification and phylogeny of the BrrHMAs

A total of 14 HMAs in turnip were identified based on the published genome results (Table 1). The HMA subfamily contains more genes in turnip (8 members) than in Arabidopsis. To determine the phylogenetic relationships between turnip and Arabidopsis HMAs, we constructed a phylogenetic tree based on HMA subfamily protein sequences. According to the orthologous relationships, the turnip HMAs were designated as BrrHMA1 to BrrHMA8 (Fig. 1). Interestingly, we found that turnip has multiple genes homologous to AtHMA2, AtHMA3, AtHMA4, AtHMA5, AtHMA6, and AtHMA7 in turnip (Fig. 1). The fourteen BrrHMAs were divided into two subclasses: the Zn/Cd/Co/Pb subclass (BrrHMA1–BrrHMA4.2) and the Cu/Ag subclass (BrrHMA5.1–BrrHMA8) (Fig. 1).

Table 1

Turnip HMA genes and their characteristics.

Gene	Gene locus	Coding region length (bp)	CDS length (bp)	Protein size (aa)	MV (kD)	pI	GRAVY
BrrHMA1	A01:730743..734744	4002	2454	817	87.63	7.65	0.141
BrrHMA2.1	A01:3562263..3566409	4147	2619	872	94.19	6.69	0.048
BrrHMA2.2	A03:27484400..27489846	5447	2724	907	97.45	7.64	−0.017
BrrHMA3.1	A03:27515852..27518987	3136	2295	764	81.98	5.82	0.252
BrrHMA3.2	A01:3551596..3554748	3153	2277	758	81.39	6.08	0.258
BrrHMA4.1	A06:15780044..15786985	6942	3573	1190	129.62	6.45	−0.268
BrrHMA4.2	A07:1087964..1095498	7535	3861	1286	138.92	7.69	−0.242
BrrHMA5.1	A09:7170712..7175590	4879	3000	999	108.55	5.85	0.144
BrrHMA5.2	A10:5078066..5081313	3248	2934	977	106.11	6.13	0.155
BrrHMA6.1	A08:12745971..12747070	1100	639	212	20.89	6.89	−0.045
BrrHMA6.2	A01:2102305..2106816	4512	2820	939	98.70	8.76	0.106
BrrHMA7.1	A06:25081334..25085063	3730	2994	997	106.60	5.08	0.254
BrrHMA7.2	A09:13230519..13235125	4607	3009	1002	107.20	4.95	0.232
BrrHMA8	A10:9935676..9941505	5830	2658	885	94.25	5.76	0.151

Fig. 1

Phylogenetic relationship of HMA proteins in Turnip and Arabidopsis. The protein sequences were aligned by ClustalW, and phylogenetic trees were constructed using the MEGA 7.0 software with the neighbor-joining method and 1000 bootstrap test. Diamonds with red and green colors represent the two subclasses of HMA proteins.

Gene structure and chromosomal location of BrrHMA genes

Turnip HMA gene lengths (including introns and exons) range from 1100 bps (BrrHMA6.1) to 7535 bps (BrrHMA4.2); their coding sequences (CDS) range from 639 to 3861 bps (Table 1), which encode 212 to 1285 amino acids (Table 1). Comparative analyses of CDS and gene sequences indicate that the BrrHMA genes contain differential exons and introns (Fig. 2). The number of exons range from 4 (BrrHMA5.2 and BrrHMA6.1) to 17 (BrrHMA6.2 and BrrHMA8) (Fig. 2). Furthermore, the number of exons do not show regular changes within two subclasses (Fig. 2).

Fig. 2

Gene structure of turnip HMA gene coding regions.

The results of chromosomal location showed that the fourteen BrrHMA genes were located on seven out of the ten turnip chromosomes (Fig. 3). Among them, chromosome A01 contained the maximum number of 4 BrrHMA genes (Fig. 3). Chromosomes A03, A06, A07, A08, A09, and A10 contained 1–2 BrrHMA members, respectively, whereas chromosomes A02, A04, and A05 carried no BrrHMA genes (Fig. 3).

Fig. 3

Distribution of HMA genes on turnip chromosomes. A01–A10 represent the chromosome number, and the ruler on the left indicates the physical map distance among genes (Mb).

Structural characteristics and physicochemical properties of BrrHMA proteins

Protein structure analysis showed that all the BrrHMA proteins except BrrHMA6.1 contain the E1-E2 ATPase and hydrolase domains whereas the members of Cu subclass have one or two HMA domains (Fig. 4). In addition, a hydrolase 3 domain was observed in the two homologous proteins BrrHMA7.1 and BrrHMA7.2 (Fig. 4). Most BrrHMA proteins except BrrHMA6.1 contain 5–8 predicted transmembrane helixes (Fig. 4). We further analyzed the diversity of conserved motif compositions in BrrHMAs using the MEME program. A total of 15 conserved motifs (motifs 1–15) were identified within the proteins (Fig. S1 in Supporting Information). The members of the same subclass generally contain similar motifs (Fig. S1 in Supporting Information). Most BrrHMA proteins contain motifs 1, 3, 4, 5, 6, 7, 8, 10, and 11 (except BrrHMA6.1), as well as motif 2 (except BrrHMA1) (Fig. S1 in Supporting Information). Each subclass includes several relatively specific motifs. The members of Zn/Cd/Co/Pb subclass have motifs 9, 12, 13, and 15, whereas the Cu/Ag subclass BrrHMAs have motif 14 (Fig. S1 in Supporting Information).

Fig. 4

Distribution of the conserved structural domains in turnip HMA proteins. The conserved domains were analyzed using the Pfam program and the diagram was drawn with TBtools software. The transmembrane helixes were manually added according to predictions by TMHMM Server 2.0.

The predicted physicochemical properties of BrrHMA proteins (i.e., MW, pI, and GRAVY values) are shown in Table 1. The MW values ranged from 20.89 (BrrHMA6.1) to 138.92 kD (BrrHMA4.2), most of which were included within 80–110 kD (Table 1). The pI values of BrrHMAs ranged from 4.95 (BrrHMA7.2) to 8.76 (BrrHMA6.2) while the GRAVY results ranged from −0.268 (BrrHMA4.1) to 0.258 (BrrHMA3.2) (Table 1).

Cd tolerance differences among turnip landraces

After Cd treatment of 16 turnip landraces, we found that the turnip landrace KTRG-B37 had the highest germination percentage index value (GPI; 0.983) and seedling survival index value (SSI; 0.5), whereas KTRG-B45a had the lowest GPI value (0.423) and a relatively low SSI value (0.012). These findings indicate that KTRG-B37 is Cd-tolerant and KTRG-B45a is Cd-sensitive.

Expression profiles of BrrHMA genes in different tissues

Tissue-specific gene expression analysis of BrrHMAs in Cd-tolerant (KTRG-B37) and Cd-sensitive (KTRG-B45a) turnip landraces revealed that BrrHMA gene expression varies in different turnip landraces. In seedlings under control conditions, expression of BrrHMA1, BrrHMA2.1, BrrHMA2.2, BrrHMA5.2, and BrrHMA6.1 was significantly higher in leaves than in roots (P < 0.05, 0.01, or 0.001). In contrast, expression of BrrHMA3.2 and BrrHMA4.2 was significantly higher in roots than in leaves (Fig. 5B). For both Cd-tolerant and Cd-sensitive turnips, expression levels of BrrHMA3.1, BrrHMA4.1, BrrHMA6.2, BrrHMA7.2, and BrrHMA8 in roots and leaves were similar (Fig. 5B).

Fig. 5

Analysis of Cd tolerance differences among different turnip landraces and gene expression levels of BrrHMA genes in different tissues. (A) Germination percentage indexes and seedling survival indexes of different turnip landraces treated by 50 mg L−1 Cd. (B) Relative expression levels of BrrHMA genes in roots and leaves of Cd-tolerant (KTRG-B37) and Cd-sensitive (KTRG-B45a) turnip landraces. Data represent means ± standard deviation (n = 3). *, **, and ***‘indicate significant differences between root and leaf for specific genes at 0.05, 0.01, and 0.001 levels, respectively.

Sequence differences of Zn/Cd/Co/Pb-subclass BrrHMA genes between Cd-tolerant and Cd-sensitive landraces

To analyze the potential roles of the BrrHMA genes in responses to Cd treatment, we cloned the CDS of Zn/Cd/Co/Pb-subclass BrrHMA genes from Cd-tolerant and Cd-sensitive turnip landraces, respectively. In Cd-tolerant and Cd-sensitive turnip landraces, all seven Zn/Cd/Co/Pb-subclass BrrHMA genes have the same coding sequences regardless of gene length (Table S4 in Supporting Information).

Expression differences of Zn/Cd/Co/Pb-subclass BrrHMA genes between Cd-tolerant and Cd-sensitive landraces under Cd treatment

Changes in gene expression in response to Cd treatment varied for Zn/Cd/Co/Pb-subclass BrrHMA genes in both the roots and leaves of Cd-tolerant and Cd-sensitive turnip landraces (Fig. 6A). Three genes (BrrHMA1, BrrHMA2.1, and BrrHMA4.1) showed specifically up-regulated expression in roots of KTRG-B37 under Cd treatment and one common up-regulated gene (BrrHMA2.2) was detected in roots of both turnip landraces (Fig. 6B and Table S5 in Supporting Information), whereas one gene (BrrHMA1) was specifically down-regulated in roots of KTRG-B45a (Fig. 6B and Table S5 in Supporting Information). Among the seven Zn-subclass BrrHMA genes, one gene (BrrHMA4.2) was up-regulated in the leaves of Cd-tolerant plants after Cd treatment; in contrast, after the same treatment in Cd-sensitive plants, four Zn-subclass BrrHMA genes (BrrHMA2.1, BrrHMA3.1, BrrHMA3.2, and BrrHMA4.2) were down-regulated in the leaves (Fig. 6B and Table S5 in Supporting Information).

Fig. 6

Gene expression results of the Zn/Cd/Co/Pb subclass members of BrrHMA genes in Cd-tolerant (KTRG-B37) and Cd-sensitive (KTRG-B45a) turnip landraces under Cd treatment. (A) Heat map summarizing the relative expression levels of BrrHMA genes in roots and leaves of Cd-tolerant (KTRG-B37) and Cd-sensitive (KTRG-B45a) turnip landraces compared to their respective control levels. Different colors correspond to the genes' log2-transformed relative expression ratios depicted in the bar at the bottom of the figure. Data represent the means (n = 3). (B) Differentially expressed BrrHMA genes in roots and leaves of Cd-tolerant (KTRG-B37) and Cd-sensitive (KTRG-B45a) turnip landraces. Red and green arrows indicate up-regulated and down-regulated genes respectively.

Discussion

Functional studies on members of plant HMA genes have been reported in several species, including A. thaliana (Eren and Arguello, 2004, Andres-Colas et al., 2006, Laurent et al., 2016), O. sativa (Lee et al., 2007, Deng et al., 2013, Huang et al., 2016), Nicotiana tabacum (Chang and Shu, 2015), Noccaea caerulescens (Lochlainn et al., 2011), Sedum alfredii (Zhang et al., 2016), Sedum plumbizincicola (Liu et al., 2017), Silene vulgaris (Li et al., 2017b), and B. napus. However, these studies have mainly focused on HMA2, HMA3, HMA4, HMA5, or HMA9. Thus, a comprehensive understanding of roles of the HMA family remains largely limited. In the present study, we identified 14 HMA subfamily members in turnip through bioinformatics analysis. These 14 turnip HMA genes exceed the 9 HMA gene found in Arabidopsis. This larger number of HMA genes in turnips may be explained by polyploidization events followed by chromosomal reduction and rearrangement, gene duplication and loss (Yin et al., 2017). While segmental duplication of HMA gene may have occurred in turnip, the duplicated homologous blocks have yet to be determined.

Our phylogenetic analysis of turnip HMA gene sequences indicates that, similar to other species (Axelsen and Palmgren, 1998, Li et al., 2015), BrrHMA genes can be divided into two subclasses: the Zn/Cd/Co/Pb subclass and the Cu/Ag subclass. In turnip, there are seven HMA genes in the ZN/Cd/Co/Pb subclass (BrrHMA1–BrrHMA4.2) and seven genes in the Cu/Ag subclass (BrrHMA5.1–BrrHMA8). This classification is supported by our analysis of predicted turnip HMA protein features. However, turnip HMA gene size and structure did not show significant similarity within subclasses. Specifically, although the gene structure is very similar between homologous gene pairs of BrrHMA2–BrrHMA4, BrrHMA5.2 and BrrHMA6.1 have fewer exons compared to BrrHMA5.1 and BrrHMA6.2, respectively. We also found that turnip HMA genes have different numbers of exons. Exon-intron increase or decrease is commonly thought to be a result of integration and realignment of gene fragments and contributes greatly to the evolution of gene families (Xu et al., 2012).

In Arabidopsis, the Cu/Ag subclass usually has three conserved domains, E1-E2 ATPase domains (PF00122), a heavy metal-associated domain (PF00403), and a haloacid dehalogenase-like hydrolase (PF00702), whereas the Zn/Cd/Co/Pb subclass lacks the heavy metal-associated domain (Li et al., 2015). Our results show that turnip HMA protein structures are highly similar to AtHMA protein structures (Williams and Mills, 2005). In contrast, turnip HMA proteins are different from P. trichocarpa HMAs (Li et al., 2015). These results suggest that the HMA subfamily has been highly evolutionarily conserved with in Brassicaceae. We found that BrrHMA6.1 only shares one heavy metal-associated domain with the other turnip HMA members. However, the BrrHMA6.1 gene sequence is highly similar to BrrHMA6.2. Furthermore, members of the Arabidopsis HMA gene family show similar differences from other members of the family, e.g., AtHMA3 only has the E1-E2 ATPase domain (Wang et al., 2014). Thus, we designated BrrHMA6.1 as a truncated HMA member. We also found that BrrHMA7.1 and BrrHMA7.2 have a supernumerary hydrolase 3 at the end of the haloacid dehalogenase-like hydrolase, which has not previously been reported in other species. However, it is unknown whether this structure has contributed to functional differentiation in these proteins. Transmembrane helixes are a typical feature of metal transporters. BrrHMA proteins are predicted to have 5–8 transmembrane helixes. BrrHMA6.1, however, is not predicted to have a transmembrane helix. Protein motif compositions also show structural similarity within BrrHMA subclasses as well as the particularity among different members, indicating that the whole family has various functions. Overall, structural analysis of all BrrHMAs indicates that these proteins have common structural features required to bind and transport HMs. Furthermore, the differential sequence structures among members of the BrrHMA subclasses determine which metal ions they transport and/or the function of the HMAs after heavy metal binding. For example, a recent study has reported that two HMA5 paralogs in S. vulgaris, SvHMA5I and SvHMA5II, can confer Cu hypertolerance by distinct mechanisms when expressed in A. thaliana (Li et al., 2017b).

Previous work has shown that different Chinese turnip landraces have distinct abilities to tolerate and accumulate Cd (Li et al., 2016). We identified a Cd-tolerant (KTRG-B37) and a Cd-sensitive (KTRG-B45a) landrace. Interestingly, although most BrrHMA genes showed similar spatial expression patterns in these two turnip landraces, some members were differentially expressed between the Cd-tolerant and Cd-sensitive landraces. To explore the underlying genetic mechanisms of differential responses to Cd in Cd-tolerant and Cd-sensitive turnip landraces, we compared sequences of the Zn/Cd/Co/Pb subclass of BrrHMA genes in these turnip landraces and changes in their expression under Cd stress. There are no differences in the CDS of HMA genes in Cd-tolerant and Cd-sensitive turnip landraces, which indicates that these genes have been highly conserved throughout the history of turnip cultivation. Turnip HMA genes, however, showed differential expression between Cd-tolerant and Cd-sensitive turnip landraces under Cd stress. Cd induced the up-regulated of four genes (BrrHMA1, BrrHMA2.1, BrrHMA2.2, and BrrHMA4.1) in the roots of the Cd-tolerant landrace, whereas only one HMA gene (BrrHMA2.2) was up-regulated in the Cd-sensitive landrace. Furthermore, BrrHMA1 was up-regulated in Cd-tolerant plants but down-regulated in Cd-sensitive plants. One HMA gene (BrrHMA4.2) was up-regulated in leaves of Cd-tolerant plants whereas four HMA genes (BrrHMA2.1, BrrHMA3.1, BrrHMA3.2, and BrrHMA4.2) were simultaneously inhibited in the leaves of Cd-sensitive plants. This difference in intraspecific gene expression has been previously reported in different cultivars of various species (Takamatsu et al., 2015, Zhou et al., 2016). These results indicate that in Cd-tolerant turnips the genes of BrrHMA Zn/Cd/Co/Pb subclass maintain or increase expression under Cd stress compared with those in Cd-sensitive turnips. Although the mechanisms underlying differences in HMA expression patterns remain unknown, this differential expression may lead to higher tolerance to Cd stress in Cd-tolerant landraces.

Conclusions

In the present study, we identified 14 candidate HMA genes from the turnip genome and analyzed the phylogenetic relationships, gene structure and chromosome distribution, as well as predicted conserved domains and motifs. Compared to the HMA gene subfamily in A. thaliana, the turnip HMA gene subfamily appears to have undergone an evolutionary expansion. According to our phylogenetic tree, the BrrHMAs are divided into two subclasses: the Zn/Cd/Co/Pb subclass and the Cu/Ag subclass. The BrrHMA members show similar structural characteristics within subclasses. To explore the potential roles of BrrHMAs in the high-Cd accumulative species turnip, we identified a Cd-tolerant landrace (KTRG-B37) and a Cd-sensitive landrace (KTRG-B45a). We subsequently analyzed the tissue expression patterns of the BrrHMA genes in both turnip landraces. Our results indicate that most BrrHMA genes show similar spatial expression patterns in the two different turnip landraces, although some members have differentially tissue-specific expression between Cd-tolerant and Cd-sensitive plants. We then conducted comparative analyses of gene sequences of the Zn/Cd/Co/Pb subclass BrrHMA genes and their expression changes under Cd stress between Cd-tolerant and Cd-sensitive landraces. Our results show that HMA genes share the same CDS but show differential expression between Cd-tolerant and Cd-sensitive plants. These findings indicate that the BrrHMA Zn/Cd/Co/Pb subclass gene that main stable expression or show increased expression under Cd stress may play a role in higher Cd tolerance in Cd-tolerant turnip landraces.

Author contributions

YP Yang conceived and designed the experiments. YS Wu and X Han analyzed the data. YS Wu, D Chen, X Li, and BQ Li performed the experiments. X Li wrote the manuscript. YP Yang and YH Yang revised the manuscript.

Conflict of interest statement

The authors declare that they have no potential conflict of interest.