Identification and expression analysis of the GDSL esterase/lipase family genes, and the characterization of SaGLIP8 in Sedum alfredii Hance under cadmium stress

Plant materials and stress treatments

Seedlings of the hyperaccumulator ecotype of S. alfredii were collected from an old Pb/Zn mining area in Quzhou City, Zhejiang Province, China. They were grown hydroponically in a growth chamber with day/night temperatures of 25 °C and a 16-h/8-h light/dark photoperiod. The seedlings were cultured in 1/2−strength Hoagland’s solution for 4 weeks. Subsequently, the roots of the experimental seedlings were immersed in 100 µM CdCl₂ as the stress treatment, while the seedlings of the control group were further cultured in 1/2-strength Hoagland’s solution. Roots, stems and leaves were sampled at 0 h, 6 h and 7 d. Three biological repeats per sample were taken at each time point and stored in a −80 °C refrigerator for subsequent use.

Identification of SaGELPs in S. alfredii Hance

The OrthoMCL algorithm (Li, Stoeckert Jr & Roos, 2003) was used to analyze 16 species, S. alfredii and 15 other related species, Phalaenopsis equestris, Oryza sativa, A. thaliana, Populus trichocarpa, Amborella trichopoda, Rhodiola crenulata, Kalanchoe fedtschenkoi, Medicago truncatula, Brassica rapa, Solanum lycopersicum, Daucus carota, Coffea canephora, Nelumbo nucifera, Macleaya cordata and Ananas comosus. The HMMER search (https://www.ebi.ac.uk/Tools/hmmer/) (Finn, Clements & Eddy, 2011) was conducted to identify and screen for possible SaGELPs containing Lipase_GDSL (Pfam: PF00657) domains in S. alfredii (Chen et al., 2018). The GELP homologous sequences in A. thaliana were obtained from TAIR (https://www.arabidopsis.org/index.jsp). All of the candidate SaGELPs were further confirmed by SMART (http://smart.embl-heidelberg.de/) (Letunic & Bork, 2018) according to AtGELPs. The basic information for SaGELPs were predicted using ExPASy (https://web.expasy.org/protparam/) (Mariethoz et al., 2018), including molecular weights, amino acid numbers and isoelectric point values.

Multiple sequence alignments and phylogenetic analyses

All of the validated SaGELP and selected AtGELP protein sequences were aligned with ClustalX in MEGA5 using GONNET as the protein weight matrix, with a gap opening penalty of 10 and gap extension penalty of 0.1. Phylogenetic trees were constructed using the Neighbor-joining method with the following parameters: text of phylogeny = bootstrap method; number of bootstrap replications = 1,000; and gaps/missing data treatment = complete deletion. iTOL (http://itol.embl.de/upload.cgi) tools were used to modify the phylogenetic trees (Letunic & Bork, 2016).

Gene structure and conserved motif predictions

Gene structure diagrams of SaGELPs were obtained from Gene Structure Display Server 2.0 (http://gsds.cbi.pku.edu.cn/) (Hu et al., 2014). Conserved motifs in SaGELPs sequences were identified using Multiple Expectation Maximization for Motif Elicitation (MEME, e value < 1e⁻¹⁰) services (http://meme-suite.org/tools/meme) (Bailey et al., 2015), with the following parameters: motif discovery = classic mode; number of repetitions = 0 or 1 occurrence per sequence; maximum number of motifs = 50; and optimum motif width = 6–100 residues. The consensus blocks in conserved domains were constructed using WebLogo (http://weblogo.berkeley.edu/logo.cgi) (Crooks et al., 2004).

Analysis of cis-regulatory elements from promoters

The cis-regulatory elements in the promoters were predicted in the 1.5 kb upstream regions of all SaGELP genes using the online website PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002). Different cis-responsive elements in the promoters were presented using RAST (http://rsat.eead.csic.es/plants/feature-map_form.cgi).

Co-expression network construction

A weighted gene co-expression network analysis (WGCNA) was used to construct co-expression regulatory networks based on profiles of differentially expressed gene responses to Cd stress, as described by Han et al. (2016). The Pearson’s correlation coefficient of the Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) value of each gene pair was calculated using the R programming language, with the correlation coefficient threshold set to 0.30 (Han et al., 2016). We screened the members of the SaGELP family and identified hub genes in the co-expression network (Langfelder, Mischel & Horvath, 2013). All eligible edges were classified according to their annotations, and we further analyze their associations with hub genes (Lotia et al., 2013). Finally, the co-expression subnetwork was visualized with Cytoscape v3.6.1 (Shannon et al., 2003).

Total RNA isolation and expression analysis

Total RNA of S. alfredii treated with 100 µM CdCl₂ was extracted from all roots, stems and leaves, using an RNA extraction kit (NORGEN, Thorold, ON, Canada). RNase-free DNaseI (New England BioLabs, Ipswich, MA, USA) was used to process genomic DNA and digest all samples. PrimeScript™ RT Master Mix (TaKaRa, Dalian, China) (Stephens, Hutchins & Dauphin, 2010) was used to produce the first-strand cDNA, which was stored at −80 °C for later use.

Quantitative Real-Time PCR (qRT-PCR) reactions were carried out using the SYBR^® Green premix Ex Taq™ (TaKaRa) reagent on the thermal circulator of an Applied Biosystems 7300 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) (Chen et al., 2018). Sequences of primers used in qRT-PCR are shown in Table S1. The relative expression level of each SaGELP gene was calculated based on the comparison threshold period (2^−ΔΔCT) method, using SaUBC9 as an endogenous reference gene (Sang et al., 2013). The heat map of the relative expression levels was constructed using online software at OmicShare (http://www.omicshare.com/). The qRT-PCR products of the expected size were analyzed by 1.5% agarose gel electrophoresis.

Heterologous expression of SaGLIP8 in yeast

The specific primers SaGLIP8-F/R were used to amplify the open reading frame of SaGLIP8 (Table S1). The purified PCR product was first inserted into the entry vector pDONR222 (Invitrogen, Carlsbad, CA, USA), and then yeast expression vector pYES-DEST52-SaGLIP8 were constructed by gateway LR reaction. The empty vector pYES2.0 was used as a control. Two vectors, expression vector pYES-DEST52-SaGLIP8 as well as empty vector pYES2.0, were transformed into the Cd sensitive mutant strain Saccharomyces cerevisiae (Δycf1) using the lithium acetate method (Liu et al., 2016). Positive colony selection was performed in the solid medium with 50 µg ml⁻¹ ampicillin and PCR reaction. The selected positive clones in the yeast liquid were cultured to an OD600 value of 0.8–1.0, and then spotted on SG-U (synthetic galactose-uracil) solid medium containing concentrations of 0, 15 and 30 µM CdCl₂. The strains in SG-U liquid medium were diluted (OD₆₀₀ = 100, 1/10, 1/100, 1/1000, 1/10000 and 1/100000), then incubated in a 28 °C incubator for 3 d (Chen et al., 2017; Liu et al., 2016). In addition, two transformed yeast cell strains were cultured on liquid SG-U medium containing 30 µM CdCl₂ for 96 h at 28 °C to determine the Cd accumulation levels by the Inductively Coupled Plasma Mass Spectrometry (ICP-MS, NexIon 300D, Perkin Elmer, Shelton, CT, USA).

Eighty SaGELP family members were identified and classified into three clades

A total of 80 SaGELPs and 56 pseudoenzymes (incomplete domain structure) were dug out (Table S2). All of the characteristics of the 80 SaGELP candidate genes are listed in Table 1, including the amino acid lengths, molecular weights and theoretical isoelectric point values. The coding sequence (CDS) lengths ranged from 900 bp (Sa28F.37) to 1,920 bp (Sa9F.272), with an average length of 1,131 bp. In total, 104 candidate sequences were obtained from A. thaliana through preliminary research, and the HMM analysis confirmed 101 AtGELP sequences and 4 pseudoenzymes (Table S3). We further used AtGELPs and SaGELPs to build the phylogenetic tree. The SaGELP gene family is divided into clades I, II and III, and the numbers of subclades are 13, 6 and 2, respectively (Fig. 1), consistent with a previous study of terrestrial plant AtGELPs (Volokita et al., 2011), showing that three branches, two major and one minor, existed in its phylogenetic tree.

Gene structure and conserved motif analysis

The genomic sequences ranged from 1,226 bp (Sa24F.223) to 3,312 bp (Sa79F.143) (Table S4). The two genes (Sa116F.65 and Sa99F.11) had maximum number of exons (seven), while the three genes (Sa39F.196, Sa9F.272 and Sa9F.81) contained only two exons. The gene structure analysis showed that the average number of exons was five, within 56 SaGELP genes (70%) containing five exons and four introns (Fig. 2B, Table S5).

Conserved sequences and motifs represent important sites for enzymatic functions. Among the 23 discovered motifs (Fig. 2C, Table S6), we analyzed four conservative motifs in the SaGELPs, blocks I, II, III and V (2, 5, 7 and 1, respectively), and the different blocks contained different motifs (Fig. 3, Table S7). A total of 13 well-conserved motifs were found (E values < 1e⁻¹⁰⁰) in most SaGELP genes, while other motifs were specific to individual SaGELPs (Table S7). Motifs 14 and 18 were only found in clade I, while motifs 16 and 21 were only discovered in clade II. Motif 23 was unique in clades I and II, and motifs 20 and 22 were distributed in clades I and III. The others motifs existed in all clades (Fig. 1A).

Analysis of cis -regulatory elements from promoters

The 1.5 kb upstream regulatory regions of the SaGELP genes were explored for stress-related regulatory elements (Table S8). The cis-acting element analysis of all genes is shown in Fig. 4. We identified cis-regulatory elements related to hormones, such as auxin, gibberellin, methyl jasmonate and ethylene. TC-rich (ATTTTCTCCA) repeats are related to cis-acting elements involving in defense and stress response. Meanwhile, some elements are also related to abiotic stress, such as heat shock elements (heat stress) and MYB-binding (AACCTAA, MRE) sites (drought-inducibility).

Co-expression network of SaGELPs

A large number of hub genes regulate potential target genes, including those related to general tolerance mechanisms and responses to Cd stress. Here, a total of 13 hub genes related to SaGELPs were obtained, as well as potential edge genes. The co-expression regulatory network involved 13 hub genes and 5 regulated different functional groups from GO (Gene Ontology) annotation (Table S7). Most of the co-expressed genes are involved in metabolic processes, growth and development, catalytic activity and biological regulation, indicating that SaGELPs have multiple functions in plants. We selected the edge genes involved in Cd tolerance from several regulatory networks, including cell wall and defense function, lipid and esterase, stress and tolerance, transport and transcription factor activity (Table S7).

As shown in Fig. 5, the major categories were transport (254 edges), transcription factor (112 edges), lipid and esterase (63 edges), cell wall and defense function (24 edges), and stress and tolerance (6 edges). The hub gene Sa0F.898 had the largest module in the Cd response gene co-expression network, with 128 nodes, including 60, 29, 19, 7 and 3 nodes related to transport, transcription factor, lipid and esterase, cell wall and defense function, and stress and tolerance, respectively. Other hub genes were also associated with different biological functions. For example, Sa13F.102 was mainly related to lipid, esterase, cell wall and defense function, while Sa26F.146 was mainly involved transport function. Therefore, SaGELPs might be involved in the induction of stress signals and function by activating transcription factors to regulate genes involved in metal transport. In addition, they was related to the enhancement of the plant’s resistance to heavy metals.

Tissue expression patterns and Cd response profiles

We used qRT-PCR to understand the functions of the SaGELP genes in S. alfredii and the tissue expression pattern under Cd-stress conditions at three time points (0 h, 6 h and 7 d). The tissues expression profiles of the genes were converted into a heat map on the basis of their expression levels (Fig. 6). All of the expression levels of SaGELP genes could be divided into the following three cases: (1) 75 genes expression significantly up-regulated at 6 h and decreased at 7 d in root (such as Sa28F.36, Sa5F.25 and Sa46F.20), stem (such as Sa27F.42, Sa45F.55 and Sa314F.6) and leaf (such as Sa13F.118, Sa0F.41 and Sa42F.134); (2) 7 up-regulated expression in stems and leaves (such as Sa10F.217); and (3) 12 down-regulated trends in roots, stems and leaves (such as Sa105F.31, Sa7F.458 and Sa184F.22).

Furthermore, the gene expression levels were greatly different in untreated samples (without Cd treatment). We performed a data analysis on the hub genes in roots, stems and leaves. Sa13F.102, Sa28F.36 and Sa5F.25 were constitutively expressed at relatively high levels in root (Fig. 7A), while the three most highly expressed SaGELPs in the stem were Sa13F.102, Sa5F.25 and Sa29F.188.1 (Fig. 7B), and the three most highly expressed SaGELPs in the leaf were Sa13F.102, Sa95F.131 and Sa29F.188.1 (Fig. 7C). Sa13F.102 had the highest expression level in all tissues. Meanwhile, the results from qRT-PCR gel image were in accordance with gene expression levels (Fig. S2).

SaGLIP8 heterologous expression enhanced Cd tolerance and accumulation in yeast

Due to the recent relationship between the Sa13F.102 gene and At5G45670.1 (AtGLIP8) from the above phylogenetic tree (Fig. S1), we designated Sa13F.102 gene as SaGLIP8. As a hub gene in the co-expression network, Sa13F.102 (SaGLIP8) was selected for functional verification for its strong induction in response to Cd stress, which implied vital roles in the Cd response in all three tissues. SaGLIP8 gene was expressed in the Cd sensitive mutant strain Saccharomyces cerevisiae (Δycf1). The SaGLIP8-overexpressive yeast grew better than the pYES2.0 yeast on a medium containing 15 and 30 µM CdCl₂, suggesting that the SaGLIP8 gene could increase Cd tolerance in yeast (Fig. 8A). Cd concentration measurements, with pYES2.0 as the control, revealed that the Cd content of SaGLIP8-overexpressive yeast was significantly greater than that of pYES2.0, and the difference between the two was extremely significant (P = 0.01) (Fig. 8B).