To identify all members of the bZIP gene family in S. plumbizincicola, HMMER3.0 was used to screen for candidate proteins in the S. plumbizincicola genome database (unpublished work) on the basis of the Hidden Markov Model profile of the bZIP domain (PF00170).¹ A BLASTP search was performed using 78 Arabidopsis protein sequences that were annotated according to previously published methods from TAIR.² Subsequently, Pfam, SMART,³ and CDD⁴ were used to confirm the presence of the bZIP domain in candidate proteins. All putative bZIP genes were named according to their homologs in Arabidopsis. The encoded protein sequences were analyzed using the online tool ProtParam⁵ to predict the amino acid composition, molecular weight, and isoelectric point (Gasteiger et al., 2005). Additionally, PSORT prediction⁶ was used to predict the subcellular localization of the proteins.

ClustalX2 was used to align the full-length SpbZIP and AtbZIP amino acid sequences. Phylogenetic trees were constructed using the maximum-likelihood criteria in MEGA 5.0, with 1,000 bootstrap replicates. The identified SpbZIP genes were divided into different groups according to the AtbZIP classification scheme. The phylogenetic tree was visualized using iTOL.⁷

The cis-acting elements in the promoter region 2 kb upstream of the SpbZIP genes were identified and then submitted to the PlantCARE database⁸ (Lescot et al., 2002). The position of the identified elements was graphically displayed using the TBtools software.⁹

The exon/intron structure of SpbZIP genes was analyzed and displayed using the GSDS platform.¹⁰ The conserved motifs in the SpbZIP proteins were identified using the MEME program (version 5.0.5),¹¹ with the following parameters: optimum motif width range of 6–50 amino acid residues and a maximum of 22 motifs (Bailey and Elkan, 1994).

The default parameters of the Multiple Collinearity Scan (MCScanX) toolkit were used to analyze gene duplication events (Wang et al., 2012). Diagrams were generated using the Circos program (version 0.69)¹² (Krzywinski et al., 2009). Non-synonymous (ka) and synonymous (ks) substitutions in each duplicated SpbZIP gene were calculated using KaKs_Calculator 2.0 (Wang et al., 2010).

Sedum plumbizincicola plants were collected from an old Pb/Zn mine in Huiping town, Quzhou city, Zhejiang province, China. The shoots from a single genotype were asexually propagated and cultivated in water in an artificial climate chamber at 25°C with a 16-h light/8-h dark cycle. The plants were grown in a half-strength Hoagland solution for about 4 weeks. Similarly growing plants were then treated with 400 μM CdCl₂. The roots, stems, and leaves were sampled at 0, 0.5, 2, 4, 8, and 12 h after the Cd stress treatment. Three biological replicates were collected for all samples.

The Total RNA Purification kit (NORGEN, Thorold, ON, Canada) was used to extract total RNA from the roots, stems, and leaves. First-strand cDNA was generated using PrimeScript™ RT Master Mix (TaKaRa, Dalian, China). The qRT-PCR analysis was performed in triplicate using the 7,300 Real-Time PCR System (Applied Biosystems, CA, United States) and the SYBR^® Premix Ex Taq™ reagent (TaKaRa, Dalian, China). Gene-specific primers were designed using the “Genes” module of the SPDE software (Xu et al., 2021a). The UBC gene was selected as the internal reference (Sang et al., 2013). Primers used are listed in Supplementary Table 5. Relative expression levels were calculated according to the 2^–ΔΔCT method (Livak and Schmittgen, 2001). The FPKM values for the SpbZIP genes were derived from the RNA-seq data (Han et al., 2016). Expression values were normalized via Z-score normalization. An expression profile heatmap was generated using the pheatmap package in R (4.0.2).

The weighted gene co-expression network analysis (WGCNA) R package was used to construct a co-expression regulatory network on the basis of the expression profiles of differentially expressed genes under Cd stress conditions (Han et al., 2016). The SpbZIP genes among the co-expressed genes with strong interconnections were designated as hub genes. The Pearson’s correlation coefficient threshold was set as 0.40 according to the FPKM values for each gene pair using the R (version 4.0.2) program (Han et al., 2016). We screened for co-expression edge genes associated with the SpbZIP hub genes and performed Gene Ontology (GO) analyses using the Gene Annotation Software for Plants (GFAP) (Xu et al., 2022). Subsequently, we classified the related genes according to their functions and visualized the relationships between nodes and edges using Cytoscape (version 3.6.1).

The SpbZIP60 coding sequence without the stop codon was fused to the mGFP-encoding sequence in the pCAMBIA1302 expression vector using the ClonExpress II One Step Cloning Kit (Vazyme, Nanjing, China). Agrobacterium tumefaciens GV3101 cells were transformed with the recombinant plasmid, which was then transferred into healthy Nicotiana benthamiana leaves for a transient gene expression analysis; the empty vector was used as a control. After co-culturing for 3 days, the leaves were soaked in a 4,6-diamidino-2-phenylindole (DAPI) staining solution to visualize nuclear DNA. The LSM 710 confocal laser-scanning microscope (Zeiss, Germany) was used to detect the fluorescence of the fusion protein.

The SpbZIP60 coding sequence was amplified by PCR and inserted into the pCAMBIA1300 vector. The recombinant plasmid was inserted into Arabidopsis (Col-0) plants via A. tumefaciens (EHA105)-mediated transformation (Zhang et al., 2006). The T₃ homozygous transgenic lines and wild-type (WT) plants were grown in a half-strength Hoagland solution. The seedlings were transferred to a solution containing 30 μM CdCl₂ after 4 weeks and grown for 7 days. The roots of the treated seedlings were immersed in a 10-mM EDTA solution for 0.5 h to remove Cd from the surface. The samples were dried and then digested with a solution comprising HNO₃ and perchloric acid (9:1 v/v) at 120–200°C in a microwave-accelerated reaction system (CEM, Matthews, NC, United States). The Cd content was determined using the 7500a inductively coupled plasma mass spectrometry system (Agilent, Santa Clara, CA, United States). Previously described 3,3’-diaminobenzidine (DAB) and nitroblue tetrazolium (NBT) staining methods were used to reveal the presence of H₂O₂ and O₂^– in situ (Chen et al., 2020). The chlorophyll content was measured according to an acetone ethanol extraction method (Li et al., 2000). Chlorophyll fluorescence was analyzed using the Dual-PAM-100 system (Walz, Effeltrich, Germany); the parameters were set, and the data were analyzed as previously described (Su et al., 2020).

Following a search of the S. plumbizincicola genome database using HMMER3.0, the identified candidate sequences were examined using CDD, Pfam, and SMART to confirm the presence of the bZIP domain (E-value < 1e^–5). A total of 92 non-redundant genes were identified as bZIP genes in the S. plumbizincicola genome. They were named according to the corresponding Arabidopsis homologs. The subsequent analysis indicated that the SpbZIP proteins comprise 117–707 amino acids (average of 303 amino acids), with a molecular weight of 13.7–77.3 kDa (average of 33.7 kDa) and a predicted isoelectric point of 5.05–10.26 (average of 7.02). Most of the identified SpbZIP proteins were predicted to localize in the nucleus, which is a characteristic of transcription factors (Supplementary Table 1).

To classify the SpbZIP genes into subgroups and elucidate the evolutionary relationships between S. plumbizincicola and Arabidopsis genes, we constructed an unrooted phylogenetic tree using the maximum-likelihood method and the protein sequences encoded by 78 AtbZIP genes and the 92 identified SpbZIP genes (Figure 1). On the basis of the phylogenetic tree, the SpbZIP genes were divided into 12 of 13 subgroups; the exception was subgroup M. There were no individual clades among the SpbZIP genes, suggesting that they were relatively conserved. Similar to the Arabidopsis homologs, most of the SpbZIP genes were classified into subgroups S and A. Subgroups J and K had the fewest genes, each with only two SpbZIP genes.

To gain insights into the structures of SpbZIP genes, their introns and exons were analyzed. Of the SpbZIP genes in subgroup S, 20 (21.7%) lacked introns. In contrast, three (3.3%) and seven (7.6%) genes contained one and two introns, respectively. Three or more introns were detected in 62 genes (68.5%) (Figure 2C). An examination using the MEME online program detected 22 conserved motifs in the SpbZIP proteins. The conserved motifs comprised 20–50 amino acids. Details regarding the 22 putative motifs are provided in Supplementary Table 2. Motif 1 (leucine zipper region of bZIP) was identified as the core conserved domain. A few subgroup-specific motifs were identified, including motifs 10 and 15 (subgroup A) and motifs 11, 12, and 16 (subgroup G). Most of the SpbZIP proteins in the same subgroup in the phylogenetic tree had common motifs, indicating a close evolutionary relationship and a high degree of conservation.

The 92 SpbZIP genes were distributed unequally among 30 S. plumbizincicola chromosomes (Figure 3). Segmental duplications of multiple genes are caused by chromosomal rearrangements (Yu et al., 2005), whereas tandem duplications, which mainly occur in the recombination region of chromosomes, usually result in the formation of a cluster of genes with similar sequences and functions (Ramamoorthy et al., 2008). During evolution, segmental and tandem duplications are the two main drivers of the expansion of plant gene families. In the S. plumbizincicola genome, eight segmental duplication events involving 16 SpbZIP genes (i.e., 17.4% of the SpbZIP genes) were detected. Among the segmentally duplicated gene pairs, SpbZIP42.1/SpbZIP42 and SpbZIP45.2/SpbZIP45.1 were distributed on chromosomes 4 and 14, respectively, whereas SpbZIP60/SpbZIP60.1 and SpbZIP61.3/SpbZIP61.2 were distributed on chromosomes 5 and 6, respectively. Additionally, SpbZIP36.1/SpbZIP36, SpbZIP53/SpbZIP53.1, SpbZIP44/SpbZIP44.1, and SpbZIP52/SpbZIP18.2 resulted from gene duplication events. Of these gene pairs, six were assigned to subgroup S. Furthermore, none of the genes were the result of tandem duplications. Thus, we speculated that segmental duplications were important for the expansion of the SpbZIP family in S. plumbizincicola. Moreover, the Ka/Ks ratios for all eight duplicated SpbZIP gene pairs were less than 0.5 (Supplementary Table 3), indicating that the SpbZIP family paralogs were primarily under purifying selection.

Next, we created two comparative syntenic maps of the association between S. plumbizincicola and Arabidopsis or Kalanchoe fedtschenkoi, which is a representative Crassulaceae plant species, to further clarify the origin and evolution of the S. plumbizincicola bZIP family (Figure 4). A total of 15 SpbZIP genes had a syntenic relationship with 17 and 48 genes in Arabidopsis and K. fedtschenkoi, respectively (Supplementary Table 4). Additionally, 20 orthologous gene pairs were detected between S. plumbizincicola and Arabidopsis, which was fewer than the 54 orthologous gene pairs between S. plumbizincicola and K. fedtschenkoi. There were more collinear gene pairs between S. plumbizincicola and K. fedtschenkoi than between S. plumbizincicola and Arabidopsis, which is in accordance with the fact S. plumbizincicola is phylogenetically closer to K. fedtschenkoi than to Arabidopsis. Some collinear gene pairs (involving 11 SpbZIP genes) among all three species were identified, implying that the orthologous gene pairs may have existed before ancestral divergence. These orthologous genes were also under intense purifying selection.

To clarify the regulatory mechanisms underlying SpbZIP expression, the cis-acting elements in the promoter sequences were analyzed using PlantCARE. The identified cis-acting elements (Figure 5) were divided into three categories (stress-responsive, plant development-related, and phytohormone responsive). The following seven abiotic stress-responsive elements were detected: ARE (important for anaerobic induction), MBS (MYB-binding site associated with drought-inducible expression), TC-rich repeat (stress-responsive element), WUN-motif (wound-responsive element), LTR (low temperature-responsive element), G-box, and W-box. At least one stress-responsive cis-acting element was detected in the promoter of all SpbZIP genes, with the exception of SpbZIP66, reflecting the importance of SpbZIP expression for plant responses to various abiotic stresses. Among the phytohormone-responsive cis-acting elements, ABRE was the most common, with 251 ABREs detected in 72 SpbZIP promoters (enrichment level of 3.49), followed by MeJA-responsive cis-acting elements (TGACG-motif and CGTCA-motif) (enrichment level of 2.49).

We used our previously published RNA-seq data to determine SpbZIP expression patterns (Han et al., 2016), which were revealed in terms of FPKM values, in the roots, stems, and leaves. The SpbZIP expression trends in the roots during the Cd treatment period were divided into four categories (Figure 6). The expression levels of 32 SpbZIP genes gradually decreased or increased over the entire treatment period. In contrast, the expression levels of 18 genes peaked at 1 day after initiating the Cd treatment. However, the genes whose expression in the roots was not induced by Cd stress had upregulated or downregulated expression levels in the stems (27) or leaves (8) in response to the Cd treatment. These results suggested that SpbZIP transcription factors may play a major role in the roots as part of the initial response to Cd stress. Transcription factors often rapidly respond to environmental cues. We further shortened and refined the treatment time and then performed qRT-PCR analysis to investigate the expression levels of 25 hub genes selected from the co-expression network. As expected, for most of the SpbZIP genes, the expression levels peaked earlier in the roots (4 h) than in the stems (8 h) and leaves (12 h) (Figure 7).

To further clarify the regulatory effects of bZIP family members on the expression of Cd-responsive genes, a co-expression regulatory network was constructed on the basis of the expression profiles of differentially expressed genes under Cd stress conditions determined in an earlier transcriptome analysis, in which 11 SpbZIP genes were annotated as hub genes. The nodes associated with hub genes were clustered according to functional categories, which reflected their association with metabolic processes, cellular activities, membranes, cells, binding, and catalytic activities (Supplementary Table 6). The Cd-responsive gene co-expression network had 189 nodes (Figure 8). The major categories included transcription factor (59 edges), transporter activity (52 edges), stimulus-response (43 edges), signaling (19 edges), and antioxidant activity (8 edges). The hub gene SpbZIP60.1 was associated with the most nodes (59), including 19 transcription factor nodes, 12 transporter activity nodes, 4 stimulus-response nodes, and 4 signaling nodes, followed by SpbZIP69.2 (34 nodes) and SpbZIP63.3 (21 nodes). Accordingly, in response to Cd stress, SpbZIP transcription factors appear to regulate the expression of downstream genes associated with diverse functions.

In this study, SpbZIP60 was one of the hub genes in the co-expression regulatory network, and its expression level was significantly upregulated in the roots during the Cd stress treatment. Hence, the subcellular localization of SpbZIP60 was analyzed to elucidate the potential functions of bZIP transcription factors in S. plumbizincicola. The control GFP signal was distributed throughout the cell, whereas the fluorescence of the SpbZIP60-mGFP fusion protein was detected only in the nucleus (Figure 9). Thus, SpbZIP60 likely functions as a nuclear protein that regulates transcription.

To further explore the function of SpbZIP60 under Cd stress conditions, transgenic Arabidopsis plants overexpressing SpbZIP60 were generated. The T₀ transgenic lines were verified by PCR using genomic DNA as the template. After analyzing the SpbZIP60 expression levels by semi-RT-PCR, the transgenic lines were cultivated to produce the homozygous T₃ lines used for the subsequent analyses (Supplementary Figure 1).

Leaf chlorosis and damages to the photosynthetic apparatus are observable symptoms of Cd toxicity. The degree of chlorosis in leaves at 7 days after initiating the Cd stress treatment was higher in the WT plants than in the SpbZIP60-overexpressing plants (Figure 10A). Histochemical staining revealed that less H₂O₂ and O₂^– accumulated in the transgenic Arabidopsis lines (OE#5 and OE#8) than in the WT control following the Cd treatment (Figures 10B,C). Meanwhile, the total chlorophyll content of the SpbZIP60-overexpressing plants where significantly higher than those of WT (Figure 10D). Chlorophyll fluorescence properties, which reflect the photochemical processes of PSII, are a useful indicator of the effects of heavy metal stress, especially Cd stress, on the photosynthetic apparatus. In the WT Arabidopsis plants, the Fv/Fm decreased, which was indicative of photoinhibition. Moreover, the inactivation or destruction of PSII resulted in an increase in the initial fluorescence (F₀). Additionally, the relative PSII electron transport rate was higher in the SpbZIP-overexpressing plants than in the WT plants (Supplementary Table 7). These results suggested that in response to Cd stress, the photosynthetic apparatus was damaged less in the SpbZIP60-overexpressing plants than in the WT plants. Next, we analyzed the Cd concentrations in hydroponically grown SpbZIP60-overexpressing lines. The Cd concentrations in the leaves and roots decreased substantially in the transgenic lines (Figure 10E). Therefore, SpbZIP60 significantly decreased the Cd concentration in the roots of the transgenic Arabidopsis plants, likely by inhibiting Cd uptake.