Genome-wide identification and transcriptional expression analysis of superoxide dismutase (SOD) family in wheat (Triticum aestivum)

Identification of wheat SOD gene family members

A computer-based methods were used to identify members of the SOD gene family from wheat reference genome IWGSC RefSeq v1.1 (https://wheat-urgi.versailles.inra.fr/Seq-Repository/Assemblies). A total of eight Arabidopsis SODs (AtSODs), twelve maize SODs (ZmSODs), and eight rice SODs (OsSODs) protein sequences were retrieved from the Arabidopsis Information Resource (TAIR10) database (http://www.arabidopsis.org/index.jsp), the Maize Genetics and Genomics Database (MaizeGDB) (https://www.maizegdb.org/), and the Rice Genome Annotation Project (RGAP) database (http://rice.plantbiology.msu.edu/), respectively. This information was used to identify SOD genes in wheat. Two methods were utilized to search the wheat protein sequences. The first method used a Hidden Markov Model (HMM) to search against wheat protein sequences and the second method used BLASTP (E-value <1e−5) to search SOD proteins compared against the wheat genome, followed by Pfam (v31.05) (http://pfam.sanger.ac.uk/search) to test whether the obtained sequences contained a SOD specific structurally conserved domain and to determine the number of SOD gene family members.

Chromosomal locations and syntenic analyses

The wheat genome GFF3 gene annotation file was obtained from the wheat database IWGSC v1.1 and the gene annotation of wheat SODs (TaSODs) was extracted from the GFF3 file. The start and end location information of TaSODs in corresponding to chromosomes was used to produce a physical map using MapInspect software.

Proteins characterization of predicted TaSODs

Characterization analysis of TaSODs were performed by using the protein identification and analysis tools on the ExPASy server10 (https://prosite.expasy.org/) (Artimo et al., 2012). The characteristics of protein length, isoelectric point (pI), molecular weight (MW), instability index, atomic composition, and amino acid composition were predicted. The online tools TMHMM (http://www.cbs.dtu.dk/services/TMHMM/) and SignalP4.1 (http://www.cbs.dtu.dk/services/SignalP/) were used to predict transmembrane domains and signal peptides of TaSODs (Nielsen, 2017). Subcellular localization prediction of TaSODs was performed using Plant-mPLoc (http://www.csbio.sjtu.edu.cn/cgi-bin/PlantmPLoc.cgi) (Chou & Shen, 2010). TaSOD members were modeled three-dimensionally using the Phyre2 (http://www.sbg.bio.ic.ac.uk/phyre2/html/) server in the intensive mode (Kelley et al., 2015).

Phylogenetic analyses of TaSODs

Phylogenetic relationship were inferred using a Maximum Likelihood (ML) method based on the LG model with MEGA7.0 software (Kumar, Stecher & Tamura, 2016). A midpoint rooted base tree was produced using the Interactive Tree of Life (IToL, version3.2.317, http://itol.embl.de) (Fang et al., 2019).

Analysis of TaSODs motifs and gene structures

Annotation information of TaSODs was interpreted using GSDS version 2.0 (http://gsds.cbi.pku.edu.cn/index.php) to produce TaSODs gene structure, intron/exon distribution, and intron/exon boundaries (Hu et al., 2014). Conserved TaSODs gene sequences were identified using a MEME suite analysis (version 4.9.1) and MAST Primer Search (http://meme-suite.org) software tools (Bailey et al., 2015). The parameters were established using known SOD protein sequences, including AtSODs, OsSODs, and ZmSODs, and were then applied to identify conserved TaSODs as follows: each sequence comprised any number of non-overlapping occurrences of each motif, the number of different motifs was 20, and motif width ranged from 6 to 50 amino acids. The functions of these predictive motifs were analyzed using InterPro (http://www.ebi.ac.uk/interpro) and SMART (http://smart.embl-heidelberg.de/), and TBtools software (https://github.com/CJ-Chen/TBtools) was used for graphical visualization.

Multiple conditional transcriptome analysis of TaSODs

Original RNA-seq data from multiple conditional transcriptome analyses were download from the NCBI data base and mapped to the wheat reference genome which reference by hisat2. After this, genes were assembled using Cufflinks to assess expression levels of TaSODs (normalized fragments per kilobase of the exon model per million mapped reads). The R package “pheatmap” was used to produce a heatmap of expression profiles for TaSODs (Song et al., 2019).

Wheat seedling growth and stress treatments

Seeds of Emai 170 (a hexaploid common wheat variety) were surface sterilized with 1% hydrogen peroxide, thoroughly rinsed with distilled water, and germinated for 2 days in a 20 ° C incubator (Ma et al., 2016). According to the reported method, the seedlings were transferred and cultured in a continuously ventilated 1/4 strength Hoagland nutrient solution (Zhu et al., 2015). After 3 days, the Hoagland solution strength was increased to 1/2. After a further 3 days, the plants were treated with 150 mM sodium chloride (NaCl), 20% polyethylene glycol (PEG) and 180 mM mannitol. The pH of the nutrient solution was adjusted to 6.0 every 2 days by using 1 M KOH or 0.2 M H₂SO₄. During the treatment, the growth environment was set at 26 ° C and 16 hours/8 h (day/night). Three biological replicates were set per treatment. Leaves and roots were collected at 2, 4, 8, 12, 24, 72, 120 h after treatment with three seedlings that were biologically mixed together. Each treatment consisted of three replicates, each of which included three technical replicates. Next, the sample was immediately frozen with liquid nitrogen and stored at −80 °C.

Real-time quantitative PCR and data analysis

To elucidate the developmental and tissue-specific expression profiles of the SOD gene in wheat, quantitative real-time PCR (qRT-PCR) was performed to detect the expression level of the SOD gene. Total RNAs from different tissues (leaf and root) and stress-treated leaves were reverse transcribed with 5X All-In-One RT Master Mix (Perfect Real Time) kit (ABM, Canada) into cDNAs for qPCR analysis. Gene-specific primers were designed using Primer Premier 5.0. For qRT-PCR assay, the cDNA was diluted to 400 ng/µL with ddH₂O. The reaction system contained 5 µL of 2 ×SYBR green Mix, 0.5 µL of each primer (10 µM), 0.5 µL of template (about 400 ng/µL) and 3.5 µL of ddH₂O to make a total volume of 10 µL. The protocol was carried out as follows: pre-denaturation at 94 °C for 3 min (step 1), denaturation at 94 °C for 10 s (step 2), primer annealing/extension and collection of fluorescence signal at 60 °C for 30 s (step 3). The next 40 cycles start at step 2. Three biological replicates were performed for each sample, with three technical replicates repeated each. Relative expression levels were determined using the 2^−ΔΔCt method (Livak & Schmittgen, 2001). The expression level of TaSOD genes were plotted using Origin software.

Identification of SODs from the wheat genome

In order to identify wheat SOD proteins (TaSODs), twenty-eight known SOD proteins, including eight AtSODs, twelve ZmSODs and eight OsSODs (Kliebenstein, Monde & Last, 1998; Dehury et al., 2013; Krishna et al., 2014), were collected as query sequences to a conduct BLASTP searches against the wheat reference genome IWGSC v1.1 (E-value < 1e–5). Candidate hits were further confirmed using Pfam and local BLASTP searches with the core motif (E-value < 1e–5) to further confirm whether the TaSODs contained the SOD domain. A rigorous bioinformatic screening identified 26 reliable TaSODs (Table 1), including eleven, five, and ten loci from the sub-genomies of A, B, and D, respectively. At this point, fifty-four SODs proteins were obtained from the four plant genomes (Arabidopsis, rice, maize, and wheat), as shown in detail in the supplemental information (Table S1). The sequences were re-named in ascending order based on their phylogenetic relationships of the OsSODs (Liu et al., 2018). The analysis of 26 wheat SOD genes produced 17 Cu/Zn-SODs (TaSOD1.1a-TaSOD1.11b), 6 Fe-SODs (TaSOD2.4-TaSOD2.9), and 3Mn-SODs (TaSOD2.1-TaSOD2.3). This result was consistent with the protein annotation information. Furthermore, alternative splicing isoforms were observed in TaSOD1.1, TaSOD1.5, TaSOD1.6, TaSOD1.7, TaSOD1.8 and TaSOD1.11.

Gene structure and chromosomal distribution of wheat genes encoding SOD proteins

In order to investigate the gene structure of TaSODs, we analyzed their GFF3-formatted annotation and observed that all TaSODs had introns. A sequence alignment of 26 TaSODs using DNAMAN software revealed low homology between the 26 proteins, and the highly conserved region was mainly concentrated at the C-terminus, which may be the key region for the function of TaSODs (Fig. 1). Exon-intron structural diversity frequently plays a key role in the evolution of gene families and can provide additional evidence to support phylogenetic grouping (Qu & Zhu, 2006; Liu, White & Macrae, 2010). The exon-intron structure of the TaSOD genes was further examined based on its evolutionary classification. As shown in Fig. 2B, all TaSOD genes contained introns, and number of introns ranged from four to sever. Sever TaSOD members (TaSOD1.9, TaSOD1.10, TaSOD1.11a, TaSOD1.11b, TaSOD2.5, TaSOD2.6, and TaSOD2.8) contained the largest number of introns (sever introns), while the smallest number (four introns) was observed only in TaSOD1.5b. As expected, the exon/intron distribution patterns were similar among SOD members within each clade of the phylogenetic tree. For example, the TaSOD2.1, TaSOD2.2, and TaSOD2.3 had the same number of exons/introns splicins and similar length.

Information regarding to TaSODs were extracted from the GFF3 reference file of the wheat genome to determine the chromosomal location of the TaSOD genes. Based on this extracted physical location (Supplemental information: Table S3), a chromosomal map of TaSOD genes was produced using the MapInspect software. The SOD gene map of the wheat genome showed SODs only on chromosomes 2, 4, and 7. The density of these loci was highest on chromosome 2 was higher with 38.46% of all SOD genes (Fig. 3).

TaSODs protein features

The amino acid sequences of 26 TaSODs proteins were submitted to the ExPASy server10 (http://www.expasy.org/tools/) online analytical system for analysis of biochemical characteristics such as isoelectric point (pI), relative molecular mass (MW) and instability index (Table 1). The results showed that the TaSODs had an average theoretical pI of 6.69 with a from 5.22 to 9.42. Protein length ranging from 141 to 392 amino acids with an average of 236 amino acids and an average molecular weight of 25.14396 kDa (range from 14.1218 kDa to 43.40094 kDa). In line with previous results, all Cu/Zn-SODs were acidic whereas Fe-SODs and Mn-SODs were basic or alkaline (Dehury et al., 2013; Zhang et al., 2016a). In the present study, most of the SOD1 were acidic in character, apart from two SOD1 enzymes (SOD1.6b and SOD1.7b). Most of the SOD2 enzymes were alkaline proteins except for two SOD2 enzymes(SOD2.7 and SOD2.9). The GRAVY analysis showed that the SOD2 subfamily contained only hydrophilic proteins, whereas the SOD1 subfamily comprised hydrophilic and six SOD1 (TaSOD1.3, TaSOD1.5a, TaSOD1.9, TaSOD1.10, TaSOD1.11a, and TaSOD1.11b) hydrophobin proteins. Interestingly, all Cu/Zn-SODs of wheat that were predicted to be localized in the cytoplasm were classified as acidic amino acids. All Fe-SODs were predicted to be localized in the chloroplasts and comprised mostly alkaline amino acids, whereas all Mn-SODs were composed with alkaline amino acids and located in mitochondria and were composed of alkaline amino acids.

Phylogenetic relationship analyses

To gain a better understanding of the evolutionary history and evolutionary relationships of the SOD gene family in wheat, a phylogenetic tree was generated using a ML method with the full-length amino acid sequences. The phylogenetic tree revealed that these SOD proteins could be classified into two major groups: SOD1 and SOD2. Moreover, we found that the SOD1 subfamily consisted of Cu/Zn-SODs and the SOD2 subfamily consisted of Fe-SODs and Mn-SODs. Based on the phylogenetic tree, we observed that SOD proteins within the same subfamily were clustered together, whereas Fe-SODs and Mn-SODs were divided into one sub-groups (Fig. 4), implying that Fe-SODs and Mn-SODs originated from a common ancestor (Alscher, Erturk & Heath, 2002). The SOD1 group consisted of 17 TaSODs (TaSOD1.1a to TaSOD1.11b), 3 from AtSODs, 6 from ZmSODs, and 5 from OsSODs. Similarly, the SOD2 proteins included 9 TaSODs (TaSOD2.1 to TaSOD2.9), 5 AtSODs, 6 ZmSODs, and 3 OsSODs. Moreover, we also could find that the dicot SODs (Arabidopsis) have more closely phylogenetic relationship related to monocot SODs (wheat, rice and maize) in each clade when compared with all plants.

Conserved motifs and clustering analyses of TaSODs

To investigate the evolutionary relationship of SODs in wheat, a phylogenetic tree was constructed by aligning the 26 TaSODs. These TaSODs clustered in two groups (TaSOD1 and TaSOD2), which was highly consistent with the type of their metal cofactors (Fig. 2A). To further examine structural diversity and predict the function of the TaSOD proteins, 20 motifs in TaSODs were identified using MEME software and were further annotated using Inter ProScan 5 software (Fig. 2C). Details of these 20 motifs are shown in the supplemental information: Table S2. Previous studies reported that the SOD gene family typically contains highly conserved domains involved in metal binding (Perry et al., 2010). The motifs 1, 2, 3, 6, 10, 11, and 17 together constitute the SOD conserved sites. The motifs 1, 2, 10, and 13 are Cu/Zn-SODs conserved domains, and the motifs 3, 6, 11, and 17 are conserved domains of Fe-SODs and Mn-SODs.

The same subfamily-associated pattern was observed regarding common motifs. All TaSODs in the SOD1 subfamily contained motifs 1 and 5, however, motif 5 is not a conserved domain of TaSOD. All members of SOD2 subfamily contained motifs 3 and 6. In addition, the Cu/Zn-SODs conserved domains (motif 1) were analyzed to understand the relationship between TaSOD1 and SODs in other species. An alignment of all Cu/Zn-SODs conserved domains of 17 TaSOD1 was termed a conserved motif 1. This result showed eight conserved amino acids (glycine, leucine, histidine, aspartic acid, serine, threonine, asparagine, and proline) in the Cu/Zn-SODs conserved domains motif. The conserved motif of Fe-SODs and Mn-SODs conserved domains site were termed as motif 3 and motif 6, respectively. Motif 3 had eight conserved amino acids (valine, proline, tyrosine, alanine, leucine, glutamic acid, serine, and histidine), and motif 6 included the conserved metal-binding domain “DVWEHAYY” of the Mn-SODs and Fe-SODs. Sequences, locations, and logos of the conserved motifs (motif 1, motif 3, and motif 6) in the TaSOD proteins are shown in Fig. 5. The data analyses supported our results. All of the identified wheat genes contained conserved domains of the sod family. In line with previous studies on different plant species, the SOD gene family in wheat contained characteristic amino acids, including a series of highly conserved active site residues that play roles in the sequence-specific binding of mental ions.

Functional categorization of wheat superoxide dismutase genes functional classification by GO annotation

Since superoxide dismutase function could scavenge intracellular superoxide anion radicals to protect cells from damage by reactive oxygen species (ROS) (You & Chan, 2015). Gene ontology (GO) annotation was performed to predict TaSOD functions (Gotz et al., 2008; Yin et al., 2018a; Yin et al., 2018b). Enrichment analysis showed that most TaSODs were annotated under GO terms ‘superoxide dismutase activity’ (GO:0004784), ‘oxidation–reduction process’(GO:0055114), and ‘removal of superoxide radicals’ (GO:0019430). Several TaSOD-I (e.g., TaSOD1.1 and TaSOD1.4) were annotated under ‘gluconeogenesis’ (GO:0006094), ‘glycolytic process’ (GO:0006096), and ‘response to cadmium ion’ (GO:0046686). In TaSOD-II subfamily, eg. TaSOD2.1 and TaSOD2.3, were annotated under ‘response to salt stress’ (GO:0009651), ‘response to zinc ion’(GO:0010043), and ‘defense response to bacterium’ (GO:0042742). GO enrichment analysis showed that multiple functions of TaSODs (Fig. 6, Supplemental information: Table S4). These results are in accordance with those of previous studies, which suggested that plant antioxidant is a complex regulating network that involves in substance and energy metabolism, and defense process (Tewari et al., 2008; Karuppanapandian et al., 2011).

Multiple conditional transcriptome analysis of TaSODs

We performed comprehensive microarray analyses to estimate the expression level of each TaSOD gene in different organs. RNA-seq data (Supplemental information: Table S5) from multiple conditional transcriptome analysis were downloaded from the NCBI database and mapped to the wheat reference genome using hisat2. After this, genes were assembled using Cufflinks software to assess the expression levels of TaSODs. The R package “pheatmap” was used to produce a heatmap of the wheat SOD genes. Previous studies showed that different pathways of SOD enzyme expression regulation patterns are unique and interact with each other (Dou et al., 2010). As shown in Fig. 7, the SOD gene family members were expressed in different tissues, and expression patterns differed between the SOD gene family members. The expression patterns were similar within subfamilies. The TaSODs can be classified into two groups: one group contains members that are widely expressed in numerous tissues, at different developmental stages, and under different treatment conditions, and the other group contains members that are not consistently expressed or highly induced only under certain conditions. Interestingly, we also found that most Fe-SODs were not highly expressed in various tissues and under different environmental abiotic stresses. Furthermore, in the salt stress environment, the expression levels of most Cu/Zn-SODs and Mn-SODs under salt stress conditions. In contrast, we clearly found that Cu/Zn-SODs and Mn-SODs were significantly up-regulated under drought and high temperature conditions. In particular, TaSOD1.1a and TaSOD1.4 showed the highest expression levels under drought and heat stress.

Expression analysis of wheat SOD genes in response to salt and drought

High salinity and drought are the major environmental stresses that frequently affect the growth and development of plants under various natural conditions (Xia et al., 2012). To simulate salt stress and drought conditions, we treated wheat seedlings with 0.15 M/L NaCl, 0.18 M/L mannitol, and 20% PEG at one-leaf stage. After treatment, the seedling growth status changed considerably. At 12 h, the leaves softened and turned yellow, and growth inhibition started. After 72 h, leaf moisture was lost, and the leaves became brittle (Fig. S2).

To understand how the TaSOD genes are involved in salt and drought stress responses, qRT-PCR was used to analyze the expression profiles of these genes under different conditions for different time in the leaves and roots. Our data show that under adverse conditions, there is a complex regulatory mechanism for the expression of TaSOD genes (Fig. 8). In leaves, the expressions of 8 genes under three treatments were variations, whereas most genes expressions were downregulated during early stage of PEG treatment. The expression of TaSOD1.7 was dynamic, increasing before 8 h, then gradually decreasing, and finally dropping to the lowest point. In root tissue, obviously, most TaSOD genes were up-regulated in response to mannitol stress. Among them, five members (TaSOD1.7, TaSOD1.9, TaSOD1.11a, TaSOD2.1 and TaSOD2.3) exhibited more than 2- to 5- fold decreases. During salt treatment, the expression levels of most TaSOD genes exhibited slight variations, whether leaf or root tissue. However, TaSOD1.7 gene showing obvious differential expression in response to NaCl stress (TaSOD1.7 was up-regulated in leaf, while it was down-regulated in root; Fig. 8).

Homology modeling of TaSODs

All 26 wheat TaSOD members were three-dimensionally modeled using the Phyre2 server in intensive mode (Fig. 9). Predicted models were based on following templates to heuristically maximize the alignment coverage, percentage identity, and confidence score for the tested sequences: template c2q2IB were used in TaSOD1.1a TaSOD1.1b and TaSOD1.4 modeling, template c1jkqD in TaSOD1.2, TaSOD1.3, TaSOD1.5a and TaSOD1.5b modeling, template d2c9val in TaSOD1.6a, TaSOD1.6b, TaSOD1.7a, TaSOD1.7b, TaSOD1.8b, TaSOD1.9 and TaSOD1.10 modeling; template d1srda in TaSOD1.11a and TaSOD1.11b models .And template c4c7uB in TaSOD2.1, TaSOD2.2, and TaSOD2.3 modeling; template c6bejA in TaSOD2.4 ,TaSOD2.5, TaSOD2.6, TaSOD2.7, and TaSOD2.9 modeling; template c1xreB in TaSOD2.8 modeling. The quality of models was validated using a Ramachandran plot analysis in which 80% of residues were within the permitted area, indicating fairly good structures of the models. However, it was apparent that in order to construct more reliable, and realistic models, more experimentally solved structures are required from SOD family proteins, particularly from plant SODs.

In the SOD1 subfamily, the secondary structures of modeled wheat proteins were primarily β-strands (26–41%) whereas α-helices occurred at only 3–14%. However, in the SOD2 subfamily, constituted the secondary structures of modeled wheat proteins were primarily α-helices primarily (45–60%), whereas β-strands occurred at only 8–12%. This is in line with the results reported in previous studies (Keerthana & Kolandaivel, 2015). Moreover, to assess similarity or divergence of generated models, structures were superimposed to calculate the percentages of structure coverage. The superimposed SOD1 subfamily models showed 69–100% structural coverage and the superimposed SOD2 subfamily models showed 51–89% structural coverage. In the SOD1 subfamily, we found that structural coverage of TaSOD1.1b and TaSOD1.4 was 100%. However, in the SOD2 subfamily, some models such as TaSOD2.5 (51%), and TaSOD2.6 (51%) showed low structural similarity but were above the twilight zone (30%). Taken together, our results suggest that SODs from each genome donor may either have been ancestrally similar to each other or originally divergent SODs could have been stabilized over a long domestication process resulting in changes in protein structures and functionality.