Genome-wide identification of ZmHMAs and association of natural variation in ZmHMA2 and ZmHMA3 with leaf cadmium accumulation in maize

Plant Materials and Phenotyping

Two experiments were conducted for analyzing Cd accumulation in leaves of 269 inbred lines. Detailed information of experimental design and plant growing conditions were described previously (Zhao et al., 2018). Briefly, for a controlled environment study, maize seedlings of 269 inbred lines were planted in plastic pots containing 14 kg clayey soil under natural sunlight in a greenhouse in 2015 and 2016, respectively. In order to simulate the extent of Cd pollution of urban soil in China (Li et al., 2014b), the middle-Cd treatment were applied as CdCl₂⋅2.5H₂O at 0.1 mmol kg⁻¹ (available Cd 18.8 mg kg⁻¹ in soil), while no CdCl₂⋅2.5H₂O was added to the soil for low-Cd treatment (available Cd of 3.28 mg kg⁻¹ in soil). Seedlings (18 d old) were cultured under low-Cd condition and middle-Cd conditions in a randomized complete block design with two replicates. After 15 days, the third and fourth leaves were harvested to analyze Cd accumulation in maize leaves. The traits for the pot experiment at seedling stage were labelled LCd15-Low (Leaf Cd concentration under low-Cd level in 2015), LCd16-Low (Leaf Cd concentration under low-Cd level in 2016), LCd15-Middle (Leaf Cd concentration under middle-Cd level in 2015), and LCd16-Middle (Leaf Cd concentration under middle-Cd level in 2016). For the field experiment, the inbred lines were planted in a highly Cd-contaminated soil (32.5 Cd mg kg⁻¹ soil) in Deyang city in 2015 and 2016. The field experiment was a randomized complete block design with two replications. After the seeds matured, five consecutive plants were chosen from the middle of each row, and the middle leaf below the tassel was harvested for measuring Cd concentration (LCd15-Field and LCd16-Field).

Harvested leaf samples were rinsed with tap water, washed three times with deionized water, and then dried at 70 °C in an oven until completely dried. The dry samples were ground to a powder using a pestle with liquid nitrogen. Subsequently, the Cd concentration was determined by using inductively coupled plasma-atomic emission spectrometry (ICP-MS) as described previously.

Identification of HMA family in maize and phylogenetic analysis

To identify the HMA genes in maize, the known amino acid sequences of the published cloned OsHMA2 and OsHMA3 (Takahashi et al., 2012a; Ueno et al., 2010) were used in the BLASTP program against the maize B73 reference genome (B73 RefGen_v3, https://www.maizegdb.org/). The protein was selected as a candidate protein only when passing the expected threshold of 1e−10, and amino acid sequence > 200 residues. Each annotated protein was examined for the existence of E1–E2 ATPase domain (PF00122) by SMART (http://smart.embl-heidelberg.de/). Multiple sequence alignment was performed by the ClustalW2.1 program with default parameters (Larkin et al., 2007). Transmembrane domains were predicted with internet-programs SOSUI (http://bp.nuap.nagoya-u.ac.jp/sosui/).

The phylogenetic analysis among maize, rice and Arabidopsis was conducted with the neighbor joining (NJ) method (1,000 bootstrap replicates) using ClustalW2.1 based on the full-length protein alignment. The phylogenetic tree was displayed in FigTree software. The exon-intron structure of ZmHMA genes was graphically displayed by the Gene Structure Display Server (Hu et al., 2014). The ZmHMA protein sequences were used to predict the conserved motifs by using the MEME Suite web server (http://meme-suite.org/) (Bailey et al., 2009) with the maximum number of motif sets at 10 and the optimum width of motifs from 5 to 300 amino acids.

Gene expression analysis

Expression patterns of ZmHMA genes in different maize tissues were analyzed using the genome-wide gene expression atlas of maize inbred B73 line that was reported previously (Stelpflug et al., 2016). Expression data under normal condition for the 6 tissues were combined from 79 distinct replicated samples (RNAseq) (Table S1). In addition, the responsiveness of each ZmHMA genes to Cd stress was analyzed by quantitative real-time PCR (qRT-PCR). Briefly, the 2-week-old plants of the B73 (a line with low Cd accumulation and a reference genome of maize) were grown in $\frac{1}{2}$ Hoagland’s nutrient solution amended with CdCl₂⋅2.5H₂O (200 µmol L⁻¹) for 0 h (used as control), 12 h, 24 h, and 48 h. Total RNA from the roots, stems and leaves was isolated using TRIZOL reagent (Invitrogen, USA) and then reverse transcription was performed with PrimeScript RT Reagent Kit (Takara, Japan). The primers used in the qRT-PCR experiments were designed by Primer 5.0 software and listed in Table S2. GADPH was used as a housekeeping gene for normalization in the present study. Subsequently, qRT-PCR was conducted using the SYBR premix Ex Taq kit (Takara) on an ABI 7500 Real-Time System (Applied Biosystems) as follows: 95 °C for 30 s; 95 °C for 5 s, 60 °C for 30 s, 40 amplification cycles, and then the melt curves were generated for verification of amplification specificity by a thermal denaturing step. The method of 2^−ΔΔCT was used to calculate the relative gene expression level between times (Schefe et al., 2006). The analysis included three biological replicates and three technical replicates for each sampling time.

ZmHMA2 and ZmHMA3 gene sequences and association analysis with cd accumulation in maize

In total, 103 maize inbred lines (stiff stalk, non-stiff stalk, and tropical or subtropical group) were randomly selected from a natural maize population of 269 inbred lines (Table S3), and used for an association analysis. To identify the sequence variations among candidate genes associated with Cd accumulation the ZmHMA2 (GRMZM2G099191) and ZmHMA3 (GRMZM2G175576), full sequences (5′-UTR, exons, intron and 3′-UTR) of these two genes were amplified from genomic DNA in 103 maize inbred lines using high-fidelity polymerase KOD FX Neo (TaKaRa). Three pairs of primers were designed using the B73 genome sequence as a reference and Primer 5.0 software (Table S2). A 150 bp overlap between each pair of primers was designed in the target regions. PCR products from 103 inbred lines were purified using E.Z.N.A and sequenced directly using the ABI 3730 sequencer. These sequences were aligned using MUSCLE 3.8, assembled using ContigExpress, and manually corrected using BioEdit 7.1 software (Edgar, 2004).

Nucleotide polymorphism in ZmHMA2 and ZmHMA3, including SNP and InDels, was identified and extracted (MAF ≥ 0.05). Allelic diversities with parameters of nucleotide diversity (π) and nucleotide polymorphism (θ) were calculated using DnaSP 6.0 software (Rozas et al., 2017). Tajima’s D statistics and Fu and Li’s statistical tests of ZmHMA2 and ZmHMA3 were selected to investigate the evolutionary pressure within the different regions. The significance of each DNA polymorphism associated with Cd leaf concentration of maize was calculated using a mixed linear model (MLM) implementing both population structure (Q) and kinship (K) in TASSEL 3.0 (Bradbury et al., 2007; Liu et al., 2015). The Q and K were calculated using genome wide SNPs as described by Zhao et al. (2018). The significance threshold of P < 0.01 was used for the candidate gene-based association analysis. Levels of linkage disequilibrium (LD) between pairs of two polymorphic loci were estimated using Haploview v4.2 (Barrett et al., 2004).

Subcellular localization of ZmHMA2 and ZmHMA3 proteins

We used ProtComp version 9 (http://linux1.softberry.com/berry.phtml) for comparing homologous proteins of known localization and pentamer distributions in the LocDB and PotLocDB databases to predict the putative protein subcellular localization. For verifying the prediction, the coding regions of ZmHMA2 and ZmHMA3 without the stop codon were amplified from B73 cDNA and cloned downstream of the 35S promoter in the pCAMBIA2300 vector that carried GFP. The primers used to amplify ZmHMA2 and ZmHMA3 were listed in Table S2. The PCR products were digested with BamHI, and directionally ligated into vector pCAMBIA2300 to construct the ZmHMA2-GFP or ZmHMA3-GFP fusion gene driven by a 35S promoter (Liu et al., 2010), respectively. In addition, equal volume suspensions of Agrobacterium strain GV3101 harboring Ti plasmids expressing either ZmHMA2-GFP and a plasma marker (PIP2A-RFP) or ZmHMA3-GFP and a plasma marker (PIP2A-RFP) were mixed (final OD600 = 0.6). These plasmids were individually expressed or co-expressed in tobacco (Nicotiana benthamiana) leaves by Agrobacterium-mediated infiltration. GFP and RFP signals were visualized using a fluorescence confocal scanning microscope (Nikon A1 i90, LSCM, Japan) at 72 h after infiltration.

Identification and properties of putative maize HMA genes

P_1B-type ATPases (HMAs) in maize were identified using the cloned OsHMA2 (Oryza sativa Japonica) and OsHMA3 (Oryza sativa Japonica) as protein queries in BLASTP searches against the B73 reference genome sequence. A total of 12 HMA genes were identified and were named ZmHMA1 to ZmHMA12 (Table 1) following the nomenclature of rice. These ZmHMA genes were unevenly distributed along the 5 out of 10 maize chromosomes. Four ZmHMAs were present on chromosome 2; three on chromosome 5; and only one gene each on chromosomes 1 and 9. The relatively high densities of ZmHMA genes were observed at the central sections of chromosomes including the centromeres and the pericentromeric regions (Table 1). Only ZmHMA6 and ZmHMA11 were noted at the bottom of chromosome 4. Gene structure analysis revealed that most ZmHMA genes displayed the complex exon-intron structure. Seven ZmHMA genes have more than nine exons. The predicted molecular weights of the 12 deduced ZmHMA proteins ranged from 42.3 (ZmHMA12) to 112.43 kDa (ZmHMA2) (Table 1). These results indicated that the HMA genes were highly conserved over the course of evolution since the pericentromeric regions of chromosomes have reduced recombination rates.

Phylogenetic and genomic structure analysis of maize HMA proteins

To investigate the phylogenetic relationship between the HMA gene family in maize, Arabidopsis and rice, the amino acid sequences of 12 HMA members from maize, 8 from Arabidopsis and 9 from Oryza sativa were used to construct a phylogenetic tree. Based on the phylogenetic tree, it clearly showed that all the maize HMA proteins were classified into the same corresponding categories in Arabidopsis and O. sativa, which include HMA I, HMA II, HMA III, HMA IV, HMA V and HMA VI clusters (Fig. 1). Functional annotation against gene ontology terms (GO; http://www.geneontology.org) shown that all ZmHMAs were involved in ATPase-coupled transmembrane transporter activity, cation transport and ATP binding (Table S4). ZmHMA1, ZmHMA1, ZmHMA3 and ZmHMA4 (HMA I, HMA II) were involved in Cadmium/Zinc ion transmembrane transport (GO:0070574, GO:0071577) for biological process categories. Eight of these ZmHMAs (HMA III–HMA VI) were listed as involved in Copper (Cu) ion binding (GO:0005507) for molecular function. In addition, each cluster contained at least one Arabidopsis and rice member in the HMA clusters (Fig. 1, Fig. 2A), which suggests that these homologous genes may have derived from multiple duplications during the evolution.

Due to the different lengths and abundant diversity of the functional domains of the HMA proteins, domain features were used for diversity analyses. As shown in Fig. S1, a total of 10 conserved motifs were identified in ZmHMA proteins by the protein domain analysis, including motifs 10 and 8 encoding the heavy-metal-associated domain (HMA domain), motifs 3, 5, 2 encoding E1–E2 ATPase domain, and motifs 4, 1 encoding the hydrolase domain. As expected, all the ZmHMA proteins possess the typical structures of the P_1B-type ATPases, including a E1–E2 ATPase domains (PF00122) (Fig. 2B). However, there were still obvious differences in the domain composition among these clusters. For instance, clusters V and VI had the largest number of motifs, while there only two motifs in cluster I. The HMA domain was absent in ZmHMA1 and ZmHMA12. ZmHMA1, ZmHMA2 and ZmHMA8 lacked the hydrolase domain (Fig. 2B). Analysis of transmembrane helices suggested that the proteins of cluster I had 6 to 9 predicted transmembrane domains (TMs) (Fig. S2, Table S5). Notably, ZmHMA3 and ZmHMA4 had exact similarity of gene structure on the same chromosome (Fig. 2C).

Expression analyses of ZmHMA genes

To investigate the expression patterns of ZmHMA genes in different tissues, an expression heatmap was constructed for 12 ZmHMA genes in the 13 different tissues of the B73 under normal conditions (Fig. 3A). The results indicated that the expression patterns of different ZmHMA genes varied greatly among tissues. Of them, the transcripts of ZmHMA2, ZmHMA6, ZmHMA7, ZmHMA11 and ZmHMA12 were highly and constitutively expressed in the various tissues. The ZmHMA3 transcript was expressed at a relatively high level in roots and nodes compared to other tissues. On the contrary, ZmHMA1 and ZmHMA4 were not expressed or exhibited at extremely low level in different tissues (Fig. 3A).

To determine which ZmHMA genes responded to Cd stress, the expression level of 12 ZmHMA genes was analyzed in B73 roots, stems, leaves under the Cd treatment using qRT-PCR. ZmHMA4 failed to be amplified in roots, leaves and stems. We found that the expression levels of ZmHMA1, ZmHMA3, ZmHMA5 and ZmHMA12 in the roots were significantly upregulated by the Cd treatment (Fig. 3B, Fig. 3D, Fig. 3E, Fig. 3L). Among these genes, ZmHMA3 was highly induced in roots and was upregulated more than 21-fold after 48 h of Cd treatment compared to normal growing conditions. In leaves, the expression levels of ZmHMA2, ZmHMA3, ZmHMA7 and ZmHMA12 were significantly induced by Cd treatment (Fig. 3B, Fig. 3C, Fig. 3G, Fig. 3L). Specifically, the remarkably higher expression of ZmHMA2 was found in leaves with a 15-fold increase after 12 h of Cd treatment. The other nine ZmHMA genes, ZmHMA1, ZmHMA2, ZmHMA5, ZmHMA6, ZmHMA7, ZmHMA8, ZmHMA9, ZmHMA10 and ZmHMA11, were markedly downregulated in the stems. Collectively, the data indicated that different ZmHMA genes exhibited variable levels of expression in different tissues and developmental stages of maize.

ZmHMA2 and ZmHMA3 gene polymorphisms and LD decay

Previous functional research reported that OsHMA1-OsHMA3 (Oryza sativa Japonica) belonged to zinc (Zn)/cobalt (Co)/cadmium (Cd)/lead (Pb) transport (Takahashi et al., 2012a). OsHMA2 plays an important role in root-to shoot translocation of Zn and Cd, and OsHMA3 is important for limiting the Cd in above-ground tissues (Takahashi et al., 2012b; Ueno et al., 2010). On the basis of phylogenetic analysis, the ZmHMA2 protein exhibited 70% identity with OsHMA2 and ZmHMA3, and 72% with OsHMA3. Comprehensive analysis highlighted ZmHMA2 and ZmHMA3, as candidate genes of particular interest.

To determine whether ZmHMA2 and ZmHMA3 genes were associated with Cd accumulation in leaves of maize, the natural variations within ZmHMA2 and ZmHMA3 genes were detected by determining the levels of nucleotide diversity. Sequence polymorphisms were detected among 103 maize inbred lines across 6,626 bp of ZmHMA2 and 3,966 bp of ZmHMA3, which covered 5′-UTR, exons, introns and 3′-UTR. SNPs and InDels at the ZmHMA2 and ZmHMA3 locus were identified (Tables S6, S7). From the putative genomic sequences in the 103 inbred lines, 205 SNPs and 62 InDels were detected from ZmHMA2, with one SNP and InDel every 23 and 103 bp, respectively (MAF ≥ 0.05) (Table S6). In addition, 139 SNPs and 24 InDels events were found in ZmHMA3, with one SNP and one InDel every 27 and 157 bp, respectively (Table S7). Nucleotide diversity analysis using sliding windows indicated that introns of ZmHMA2 had the highest nucleotide diversity (π = 17.58 × 10⁻³) and its exon had the lowest nucleotide diversity (π = 6.24 × 10⁻³) in maize (Table 2). However, in ZmHMA3, the estimates of nucleotide diversity in the 3′-UTR (π = 33.57 × 10⁻³) was obviously higher than that of other regions. Fu and Li’s F and D test of each region showed that the exon of ZmHMA2 and the 3′ UTR of ZmHMA3 were under purifying selection with significantly positive values. These results suggested that purifying selection and/or population size expansion had occurred in the two regions respectively. In addition, the combined LD analysis of ZmHMA2 and ZmHMA3 showed LD decay was close to 0.1 within 1300 bp (Fig. 4).

ZmHMA2 and ZmHMA3 genes associated with Cd accumulation in leaves of maize

The association analysis between natural variations in ZmHMA2 and ZmHMA3 and Cd concentration in maize were performed using the Q+K model. One co-located InDel in the intron (InDel S1174, P = 7.96 × 10⁻³) in ZmHMA2 was associated with leaf Cd concentration at seedling stage under low Cd condition in 2015 and 2016 (LCd15-Low and LCd16-Low) (Table 3). For ZmHMA3, 20 significant loci (15 SNPs and five InDels) were associated with leaf Cd concentration under various soil environments at seedling and maturing stages in 2015 and 2016 (P < 0.01) (Fig. 5, Table 3). R² marker values ranged from 7.7 to 16.3% among all those associated loci. On average, each site explained 11.2% (±3.3%) of phenotypic variance for its associated trait. Among the 20 loci, S917 located on the second exon region were significantly associated with LCd16-Middle (P < 0.01) (Table 3). Moreover, eight loci (S3097, S3114, S3297, InDel S3517, InDel S3594, S3692, S3695, S3792) with a complete LD were significantly associated with LCd15-Field and LCd16-Field. Furthermore, five SNPs and two InDels were located on the CDS region and caused the amino acid substitutions and insertions or deletions.

Putative ZmHMA2 and ZmHMA3 effect on Cd accumulation

Based on the average phenotype of plants grown in 2015 and 2016, the effects of gene polymorphisms on Cd accumulation were further analyzed using mixed linear model. Genotype with a deletion of GCA (InDel S1174) in an intron of ZmHMA2 exhibited a significantly higher Cd concentration than that of the genotype with an insertion of GCA (InDel S1174) under low Cd condition at seedling stage (Fig. 6A). The individuals carrying the minor frequency alleles had a 1.03 mg lower leaf Cd content than those with major frequency alleles (P = 8.48 ×10⁻³). In addition, the site S917 (CC) associated with LCd16-Middle at seedling stage had one amino acid substitution from alanine (A) to glycine (G) (P < 2.09 × 10⁻³). The Cd concentration of leaves in the genotype at seedling stage with SNP G/G was about 6.9 mg higher than those with SNP CC (Fig. 6B). In the field study across 2015 and 2016, four polymorphic loci (S3114, S3297, InDel S3517, InDel S3594) in the CDS region of ZmHMA3 were significantly co-associated with LCd16-Field and LCd15-Field. From the two SNP positions, two amino acid substitutions, from glycine (G) to alanine (A) at residue 3114 and from valine (V) to alanine (A) at residue 3296, were identified (Table 3). More importantly, InDel S3517 (GTCCCGTGA) encoded valine (V), proline (P) and the stop codon. In addition, InDel S3517 and the three loci above showed a complete LD pattern. A significant difference in Cd accumulation of leaves among the accessions with homozygous allelic SNP S3114, S3297 and InDel S3517 in ZmHMA3 was found at mature stage in the field study (Fig. 6C). The inbred lines with major frequency alleles (GG +TT+9) exhibited a significantly lower Cd accumulation in leaves than that of inbred lines with minor frequency alleles (CC+CC+0) (P <  2.82 × 10⁻³). The Cd accumulation in mature leaves of individuals with minor frequency alleles was significantly reduced by 10.96 mg compared to those with major frequency alleles.

Subcellular Localization of ZmHMA2 and ZmHMA3 Proteins

Possible subcellular locations of ZmHMA2 and ZmHMA3 proteins were predicted using ProtComp 9.0. This prediction suggested that ZmHMA2 and ZmHMA3 might be localized on PM (Table S8). To examine the predictions above, the control vector pCAMBIA2300-35S::GFP, the loaded vector pCAMBIA2300-35S::ZmHMA2-GFP and pCAMBIA2300-35S::ZmHMA3-GFP, were transferred into tobacco using transient transformation. ZmHMA2 and ZmHMA3 proteins were successfully expressed as fluorescent protein fusions (ZmHMA2-GFP and ZmHMA3-GFP) (Fig. 7). Compared with the control (Fig. 7D), the green fluorescence signals of ZmHMA2-GFP (Fig. 7E) and ZmHMA3-GFP (Fig. 7F) were preferentially detected in the cytoplasm of the tobacco leaf epidermis. To further confirm subcellular localization of proteins, we co-infiltrated the tobacco leaves with Agrobacterium containing the PM marker (PIP2A-RFP) with ZmHMA2-GFP (Fig. 7H) and ZmHMA3-GFP (Fig. 7I), respectively. Confocal microscopy analysis of co-infiltrated cells indicated that GFP green signals of ZmHMA2-GFP (Fig. 7K) and ZmHMA3-GFP (Fig. 7L) overlapped with the RFP signals of the PM marker. These results demonstrated that ZmHMA3 and ZmHMA3 proteins were primarily localized on the PM.