Genome-wide identification and characterization of the soybean SOD family during alkaline stress

Identification of SOD family genes in soybean genome

To identify all potential family members of SOD in soybean, the known soybean SOD amino acid sequences were used as queries to establish a Hidden Markov model (Gopavajhula et al., 2013), and searched in soybean genome database by using the HMM profile (build 2.3.2) (Finn, Clements & Eddy, 2011). The Pfam and SMART database were used to remove incomplete domains and overlapping genes (Finn et al., 2016). The molecular weight and isoelectric point values of SOD proteins were predicted using online software ExPASy (http://au.expasy.org/tools/pi_tool.html) (Artimo et al., 2012).

Bioinformatics analysis of SOD family genes

The neighbor-joining phylogenetic tree was generated by using software MEGA 5.0 (Kumar et al., 2008). The exon-intron structures were analyzed by the GSDS online server (http://gsds.cbi.pku.edu.cn/) (DuanMu et al., 2015). The conserved motifs were analyzed by MEME online server (http://meme-suite.org/) (Bailey et al., 2009). The multiple sequence alignments were performed using Clustal X program (Larkin et al., 2007). The syntenic analysis was conducted using Circos-0.69 (http://circos.ca/) (Krzywinski et al., 2009). The synonymous (Ks) and nonsynonymous (Ka) substitution rates of syntenic paralogs were calculated using DnaSP software (version 5.10.01) (Krzywinski et al., 2009). The PLANT CARE online software was used to analyses the cis-acting elements of promoters (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002).

Expression patterns of SOD family genes based on transcriptome sequencing data

To examine the expression profiles of soybean SOD family genes in different tissues, the expression profile data (GSE29163) was obtained from NCBI. The transcriptome data of soybean was downloaded from the NCBI database under cold (GSE117686), salt (GSE57252) and water-deficit (GSE49537) stresses. As there is limited information exist on soybean in response to alkaline stress, we downloaded wild soybean transcriptome data under alkaline stress. The hierarchical clustering trees of SOD family genes were generated using TM4: MeV4.9 software (Saeed et al., 2006).

Plant material, growth condition and alkaline stress treatment

The soybean (DN50) seeds were treated with 75% ethanol for 1 min and washed with sterile water before germination 2 days. Then, seedlings were grown in Hoagland nutrient solutions at 70–80% relative humidity, 22–28 °C room temperature and 8 h dark/16 h light. Twelve days after sowing, soybean seedlings were transferred into Hoagland solutions with 50 mm NaHCO₃ for 0, 6 and 12 h (DuanMu et al., 2015). The roots and leaves were harvested as three biological replicates.

Transcript data analysis of SOD family genes under alkaline stress

Total RNA was extracted from soybean using the plant total RNA isolation kit (TIANGEN, Beijing, China), and the cDNAs were synthesized using the TransScript All-in-One First-Strand cDNA synthesis SuperMix for qPCR kit (TransGen Biotech, Beijing, China). qRT-PCR assays were performed using UtraSYBR Mixture (low ROX) and ABI 7500 sequencer. The GmGADPH was used as internal control (Huis, Hawkins & Neutelings, 2010). The SOD genes and GmGADPH primers are listed in supplementary Table S1. The three biological replicates were obtained and expression levels were calculated using 2^−ΔΔCt method and Student’s t-test (Livak & Schmittgen, 2001).

Identification and phylogenetic analysis of SOD genes in soybean

In soybean, ten SOD isoenzymes were identified in different tissues (Wang et al., 2016). For example, one MnSOD isoenzyme was mainly detected in stems and seeds. Four Cu/ZnSOD isoenzymes were detected in roots, leaves, stems and seeds. While, only five FeSOD isoenzymes was mainly detected in leaves. To further explore the SOD genes in G. max, we used SOD amino acid as query sequences to search the G. max database from NCBI. A total of 19 SOD candidate proteins were obtained. All candidate proteins were subjected to Pfam and SMART database to remove incomplete domains and overlapping genes. As a result, we obtained 13 non-redundant SOD genes in soybean. To confirm the classification and evolutionary relationships of soybean SOD family members, the full-length protein sequences of GmSODs were used to construct a neighbor-joining phylogenetic tree with Arabidopsis, Medicago truncatula and Solanum lycopersicum SOD family members. The results revealed that the soybean SOD family members were classified into three subfamilies, including MnSODs (GmMSD1–2), FeSODs (GmFSD1–5) and Cu/ZnSODs (GmCSD1–6), which is consistent with previous studies (Fig. 1). In addition, compared with previous study (Wang et al., 2016), we identified three more SOD genes in soybean.

As shown in Table 1, these GmSOD family members were chosen for further protein information analyses, including the exon numbers, length of CDS, protein sequence lengths, molecular weights (MW), and theoretical isoelectric points (pI) values. The protein sequence length ranged from 152 (GmCSD1 and GmCSD6) to 314 (GmFSD1) amino acids (aa). The MW varied from 15.19 (GmCSD6) to 35.86 (GmFSD1) kDa and the pI values ranged from 5.22 (GmMSD1) to 8.57 (GmFSD2). The Cu/ZnSODs subfamily members had a lower protein sequence length and MW than others. The MnSODs subfamily members appeared with higher pI values. This finding indicated that three subfamilies may displayed a great diversity in soybean.

Phylogenetic and gene structure analysis of SOD genes

To further investigate the evolutionary relationship of SOD genes in soybean, we constructed a neighbor-joining phylogenetic rootless tree with the full-length protein sequences of GmSODs. The results suggested that three SOD subfamilies have a distant evolutionary relationship (Fig. 2A). However, we found that each SOD subfamilies have high bootstrap support pairs, such as GmFSD3 and GmFSD5, GmMSD1 and GmMSD2, GmCSD2 and GmCSD3. This indicated that each soybean SOD subfamilies was evolutionarily conserved, even though three subfamilies displayed a distant evolutionary relationship.

The diversity of gene structure is mainly influenced by the evolution of multigene families (Mercereau-Puijalon, Barale & Bischoff, 2002; Pellicer et al., 2018). To further explore the structural diversity of SOD genes, the characteristics of exon-intron structures were analyzed by the GSDS online server. As shown in Fig. 2B, different subfamilies displayed variation in exon-intron structures. For example, the MnSODs subfamily comprised of six exons as well as exhibited a similarity in genomic structure. The FeSODs subfamily members consisted of nine exons. We also noticed that GmFSD2 exhibited a longer genomic structure than 10 KB. However, three Cu/ZnSODs subfamily members (GmCSD1, GmCSD5 and GmCSD6) possessed seven exons, while others appeared with eight exons. These results further confirmed that the different SOD subfamilies diverged greatly in soybean.

Conserved domain analysis of SOD proteins

For identification of the conserved domains of SOD proteins in soybean, the conserved motifs were analyzed by the MEME online server. Our results showed that eight predicted conserved motifs were identified (Fig. 3). Among them, all proteins in Cu/ZnSODs subfamilies contained motifs 4, 7 and 8. Motifs 1, 2, 3, 4 and 5 were detected in MnSODs subfamily members. Motifs 1, 2, 3, 4, 5 and 6 were discovered in all FeSODs subfamily members, except GmFSD2. It is noteworthy that motif 1 was detected in the N-terminal domain (alpha-hairpin domain) of MnSODs and FeSODs subfamily members. Besides, motifs 2, 3, 4 and 5 were also detected in the C-terminal domain. The motifs 4, 7 and 8 were also detected in the copper/zinc SOD domain of Cu/ZnSODs subfamily members (Figs. S1 and S2). These results indicated that the MnSODs and FeSODs subfamily members showed high similarities in conserved sequences. However, the Cu/ZnSODs subfamilies displayed wide divergence with other two subfamilies, which further verify the potential diversity of functional divergence in soybean SOD genes from different subfamilies.

Chromosomal locations and syntenic analysis

To determine the genomic distribution of SOD genes, we identified their gene locations on the chromosomes and potential genome duplication events using the syntenic analysis. The thirteen SOD genes were randomly distributed among ten chromosomes, and each chromosome had one or two genes (Fig. 4). Gene duplication plays important role in plant functional diversity and evolutionary mechanism (Bowers et al., 2003). To detect the potential genome duplication events, a total of eight pairs of SOD syntenic paralogs were found in soybean genome, except GmCSD5. The results indicated that the soybean SOD family existed a high gene family expansion. Meanwhile, we found that the duplication pairs also exhibited a high evolutionary relationship in phylogenetic analysis (Fig. 2A).

The non-synonymous/synonymous substitution (Ka/Ks) rates were used to detect the positive or negative selection of gene family expansion (Kaehler, Yap & Huttley, 2017) . To further investigate whether selective pressure was associated with the SOD family genes, the Ka/Ks rates were calculated for the full-length CDS sequences of eight duplication pairs (Table S2). Previous studies reported that Ka/Ks < 1 indicates a negative selection and Ka/Ks > 1 indicates a positive selection (Wang et al., 2005). The results showed that the Ka/Ks rates for six SOD duplication pairs were more than 1. However, two SOD duplication pairs were less than 1. This finding indicated that the soybean SOD genes had experienced positive and negative selection pressure after the duplication events.

Analysis of cis-acting elements of SOD gene promoters

The cis-acting elements play significant roles in determine the regulatory roles under various stresses (Kimotho, Baillo & Zhang, 2019). In this study, to explore the potential roles of the soybean SOD genes under environmental stresses, we analyzed the SOD gene promoters coupled with 3 KB genomic sequence upstream regions of the ATG using PlantCARE online tool (Sun et al., 2014). The SOD family genes possessed a variety of putative abiotic stress and hormone-related responsive elements, including ABA (ABRE), Methyl jasmonate (MeJA), Salicylic acid, Gibberellin, MBS (Drought) and Defense and Stress responsive elements (Fig. 5; Table S3). Among them, almost all SOD family genes contained ABA-responsive elements, except GmMSD2. Ten of SOD family genes had MeJA-responsive elements. However, the GmMSD2 and GmFSD5 only contained two responsive elements, respectively. We also found that only Cu/ZnSODs subfamily genes contained Defense and Stress responsive elements. Collectively, the results suggested that SOD family genes might be related to the responses of abiotic stresses and hormone stimuli. Also, these genes might have potential functional diversity because they contained different responsive elements.

Expression profiles of SOD genes in soybean tissues

To investigate the potential function of SOD genes in soybean growth and development, we analyzed their expression in nine tissues using RNA-seq data (Chen et al., 2012). The results showed that the numbers of expressed soybean SOD genes in different tissues exhibited higher variations (Fig. 6). For example, GmMSD1, GmMSD2, GmFSD3, GmFSD5 and GmCSD1 showed higher expressions in early-maturation stage whole seed. GmCSD1, GmCSD2, GmCSD3 and GmCSD6 had higher expression levels in leaf. In addition, the GmCSD1 displayed the higher expression level in five tissues, which is consistent with the cucumber CSD1 gene studies (Zhou et al., 2017). The GmFSD2 showed an extremely lower expression levels than others in nine tissues. However, almost all the soybean SOD genes had the lowest expression level in mid- and late-maturation stage whole seed. In conclusion, these results indicated that soybean SOD genes showed special tissue expression profiles, indicating their potential divergent functions in soybean growth and development.

Expression analysis of soybean SOD genes in response to abiotic stresses

Previous studies have showed that the SOD genes play significant roles in plant responses to various stresses (Jiang et al., 2019). To investigate the potential roles of soybean SOD genes in various abiotic stress responses, we identified their expression patterns under alkaline, salt and cold stresses using transcriptome sequencing data. However, there is limited information exist on soybean in response to alkaline stress, we downloaded wild soybean transcriptome data under alkaline stress. As shown in Fig. 7, the heat map revealed that four and six of soybean SOD genes were differentially expressed (|Log₂ fold change| > 1, P < 0.05) under alkaline and salt stresses, respectively. However, only GmFSD2 displayed slightly up-regulation under cold stress. Among them, we found that GmFSD3, GmFSD5 and GmCSD5 were all up-regulated under alkaline and salt stresses, indicating that they might act as positive regulators. It is worth noting that GmFSD3 and GmFSD5 showed differential expression in soybean leaves and roots under salt stress. This finding implied that these two genes may be involved in different signaling pathways in leaves and roots under salt stress.

Expression analysis of SOD genes in response to alkaline treatment

Saline-alkali soil is considered as a major threat to crop growth and yields. However, only a few studies have focused on the mechanisms of plants respond to alkaline stress. In this study, to further confirm the soybean SOD genes in response to alkaline treatment, we detected their expression patterns in roots and leaves under 50 mm NaHCO₃ treatment by using qRT-PCR. As shown in Fig. 8A total of six genes were up-regulated expressed and only GmFSD5 gene was down-regulated in soybean leaves under alkaline treatment. Three genes (GmFSD3, GmFSD4 and GmFSD5) were significantly induced in soybean roots under alkaline treatment, whereas others showed down-regulated expression, especially at 12 h. Among them, GmFSD5 had contrary expression pattern in leaves and roots, indicating this gene may be involved in different pathways.

In addition, we found that FeSODs and Cu/ZnSODs subfamily genes displayed different roles in roots and leaves after alkaline treatment. For example, only three FeSODs subfamily genes (GmFSD3, GmFSD4 and GmFSD5) were significantly induced (P < 0.05 or P < 0.01) in roots. Five Cu/ZnSODs subfamily genes (GmCSD1/2/3/4/5) were up-regulated in leaves. However, almost all FeSODs subfamily genes didn’t express significantly in leaves, except GmFSD5 which was down-regulated expressed. To conclude, qRT-PCR confirmed the results that soybean SOD genes possibly participate in responses to alkaline.