Cloning and comparative modeling identifies a highly stress tolerant Cu/Zn cytosolic super oxide dismutase 2 from a drought tolerant maize inbred line

RNA isolation

Leaf tissue was sampled from 21-day old maize seedlings (Zea mays, HKI 335 inbred line) growing in pot under normal controlled conditions in greenhouse. and RNA was isolated using RNA isolation kit (Genetix Life Sciences, New Delhi, India). Quality and concentration of RNA was checked on Nanodrop Spectrophotometer (Thermo Fisher, Waltham, MA, USA). 1 µg of RNA was used for cDNA synthesis using oligodT primers (cDNA synthesis kit; Invitrogen, Waltham, MA, USA).

PCR amplification

Gene specific primers were designed and synthesized for maize (Sod2) (sod2_Forward ACAATGGTGAAGGCAGTTGCTGTC and sod2_Reverse TTAGCCTTGGAGCCCAATGATACC) and amplification was done by polymerase chain reaction (PCR) using high fidelity polymerase (Pfu). The amplified product was ligated into pJET vector and transformed into competent E. coli cells (DH5α strain). Colony PCR was done using above primers to identify transformed colonies and plasmid was isolated. This plasmid was sent for sequencing. BLAST analysis of sequencing results confirmed our sequences as SOD2 mRNA which were then submitted in NCBI (ALF00121.1). The sequences were retrieved for further bioinformatic analysis.

Sequence search

The target sequence was searched for similar sequence using the BLAST (Basic Local Alignment Search Tool) (https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE=Proteins) against Protein Database (PDB) (Johnson et al., 2008). Only full-length protein sequences were considered for in silico analysis.

Multiple sequence alignment and phylogenetic analysis

Similar sequences of SOD2 were searched using BLAST and multiple sequence alignment was done using Clustal W (https://www.genome.jp/tools-bin/clustalw), in order to infer the phylogenetic relationship between members of the antioxidant gene family (Thompson, Gibson & Higgins, 2003). The final alignment was completed through manual editing. The FASTA multiple sequence alignment was used to infer the neighbor-joining (NJ) phylogenetic tree using the MEGA 6.0 software (Tamura et al., 2013). The bootstrap interior branch test, superior to bootstrap analysis in testing confidence in a given tree node, was applied to test the statistical significance of the SOD2 of NJ phylogenetic tree.

Secondary structure prediction

The target sequence in FASTA format have been used for prediction of secondary structure of protein. With the utilization of online bioinformatics tool known as SABLE (http://sable.cchmc.org/), one can identify the secondary structure of target protein sequence (Adamczak, Porollo & Meller, 2005). This tool generally predicts real valued relative Solvent AccessiBiLitiEs of amino acid residues in protein. It uses evolutionary profiles for improved prediction of secondary structure. The result of secondary structure was visualized using POLYVIEW-2D.

Modelling of protein’s three-dimensional (3D) structure

The 3D structure of the protein was identified using online bioinformatics tool called MODWEB (https://modbase.compbio.ucsf.edu/modweb/) (Eswar et al., 2003). It modeled the structure using comparative protein structure modeling, which can calculate comparative models for a large number of protein sequences, using many different template structures and sequence-structure alignments.

Model validation

The model was evaluated on the basis of geometrical and stereo-chemical constraints using PDBsum (http://www.ebi.ac.uk/pdbsum) (Laskowski et al., 2018) and ProSA-Web (https://prosa.services.came.sbg.ac.at/) (Wiederstein & Sippl, 2007; Gill & Tuteja, 2010). PDBsum is used for validation of PDB structure through which stereochemical quality and accuracy of the predicted model is evaluated, which is based on Ramachandran plot calculation. Amino acid residues in proteins were examined by torsion angles Ø and Ψ and a percentage quality measurement of the protein structure was used, in which four sorts of occupancy were called core, allowed, generously allowed and disallowed regions. SuperPose (http://superpose.wishartlab.com/) was used for protein superposition that calculates the superpositions using a modified quaternion approach (Maiti et al., 2004). We calculated the root mean squared deviation (RMSD) to identify the distance between two objects by superposing the two structures and is performed between equivalent atom pairs.

Physiochemical properties

The physiochemical properties of the target protein such as molecular weight, number of amino acids, theoretical isoelectric point (pI), total number of negatively charged (Asp + Glu) and positively charged (Arg + Lys) residues, amino acid and atomic composition, extinction coefficient, estimated half-life, grand average of hydropathicity (GRAVY), instability and aliphatic index are identified using online tool ExPASy’s ProtParam (http://expasy.org/tools/protparam.html) (Gasteiger et al., 2005).

PTM site prediction

For the prediction of PTM sites, the online bioinformatics tool PTM-ssMP (Liu et al., 2018) was explored.

Protein-protein interaction analysis

For predicting the function of query sequence with the reference sequence, protein-protein interaction (PPI) networks analysis has been done using online tool STRING (http://string-db.org/) (Szklarczyk et al., 2019).

Domain and motif identification

We identified the domain composition of target protein sequence by using the bioinformatics tool SMART (http://smart.embl-heidelberg.de/) (Letunic, Doerks & Bork, 2012). It annotates and identifies the genetically mobile domains and analyses the domain structure. It can also detect extracellular, signalling and chromatin associated proteins bearing more than 500 domain families. ScanProsite server (https://prosite.expasy.org/scanprosite/) have been used for identification of signature motif in the query protein sequence (De Castro et al., 2006).

Cloning of Sod2 from HKI 335 maize inbred line

The ZmSod2 gene was amplified from cDNA using gene specific primers and cloned in standard plasmid vector, followed by sequencing the ZmSod2 HKI 335 sequence was deposited in the GenBank database of National Centre for Biotechnology Information (NCBI), USA with accession number ALF00121.1. The cloned sequence (151 amino acid) contained all the characteristic conserved domains of Cu-Zn Superoxide Dismutase superfamily (PLN02386 family protein) (Fig. S1).

Phylogenetic analysis

The multiple sequence alignment shows aligning of SOD2 from different stress tolerant plant species. Similar sequences of SOD2 from different species were aligned and further used for phylogenetic analysis using the neighbour joining (NJ) method. The results show that SOD2 in different plant species demarcated into two prominent clusters i.e., monocots (6) and dicots (5) (Fig. 1).

Secondary structure and comparative modeling

Using an online bioinformatics tool known as SABLE, we predicted the secondary structure of the target protein sequence (Fig. 2). The secondary structure of the protein sequence viewed by polyview-2D show the absence of alpha-helices. At the same time, there are a large number of coils and beta-sheets with fully exposed relative solvent accessibility (RSA). The result indicate that it is a soluble protein and no transmembrane domains have been found (Fig. S2). Most of the beta-sheets have RSA value greater than 2, which shows that the amino acids in the protein secondary structure are completely buried. The high value of the RSA of the target protein suggests that the protein is stable for a number of biological functions within the cell. The theoretical structure of SOD2 was generated using MOD web by comparative modeling of protein structure prediction through which identical or non-identical information about the target sequences is analyzed. We developed 14 models calculated for SOD2 out of which two models were selected. The model with LONGEST_DOPE had a score of −1.834 by using Modweb which was dynamically refined and validated. Comparative modeling revealed that the cloned SOD2 exhibits 83% identity with Potentilla atrosanguinea SOD2 (PDB ID: 2Q2L_A) with an e-value of 4e−86 (Figs. 3A and 3B).

Potentilla atrosanguinea, commonly known as Himalayan cinquefoil or ruby cinquefoil, is a vigorous herbaceous perennial of the rose family native to mountain slopes at lower elevations in the Himalayas. Surprisingly, SOD2 from Potentilla atrosanguinea have been previously reported as a highly thermostable SOD that can tolerate autoclaving as well as sub-zero temperatures (Kumar, Sahoo & Ahuja, 2002; US Patent US20070269811A1). Further, a single C95A amino acid substitution in the wild type Potentilla atrosanguinea SOD2 is known to enhance thermal properties (Kumar et al., 2012). These findings prompted us to further investigate the SOD2 from HKI 335 maize inbred.

Model validation

With the PDBsum, it has been shown that generated model of SOD2 protein revealed (~98.0% expected): 193 (84.6%) residues falling in most favoured region, (~2.0% expected): 34 (14.9%) residues in additionally allowed region, and 0.4% residues in generously allowed region with no residues in the disallowed region of the Ramachandran plot (Fig. 3C). Z-score of PROSA energy indicating overall model quality was used to check 3D models of protein structures for potential errors. In these plots displaying Z-scores, value (−6.51) of the target model was determined by X ray crystallography (represented in light blue) and nuclear magnetic resonance (represented in dark blue). This value was extremely close to the value of template 2Q2L (−6.44) (Fig. 3D). Root mean squared deviation (RMSD) value indicates the degree to which the two three dimensional structures are similar. RSMD analysis of the SOD2 model was measured from its template (2Q2l_A). The Cα RMSD and backbone RSMD deviation for both the target and the template was 0.27 Å (Fig. 3E).

Physiochemical properties

The target protein consists of very high percentages of glycine (18.5%), alanine (7.9%) and valine (9.9%) in comparison to other amino acids (Table 1). The high percentage of glycine (18.5%) shows that triple helical structure of the protein is likely to be more stable. Like that of glycine residues in protein, the proline also plays a major role in stability of helix of protein’s secondary structure. Here, the residues of proline show less percentage (5.3%). On the other hand, the percentage of glycine, valine and leucine among the hydrophobic groups are 32.56% and 17.44% respectively, while serine and threonine (both hydrophilic group) are 20.51% and 28.21% respectively. The isoelectric point (pI) value was 5.43, while the instability index was 20.1 (Table 2). The aliphatic index was 80.66, which is extremely high. The value of grand average of hydropathy (GRAVY) was −0.132 that lies between −2 and +2, suggesting that the protein is hydrophobic in nature and rated positively. The atomic composition of the protein is shown in Fig. 4.

Post translational modification (PTM) site prediction

The post translational modification (PTM) is described as amino acid modification on the basis of protein sequence and also considered as vital issues for regulating the physiological and biological functions inside the cell. We could identify only one PTM site (phosphorylation), located at position 20 having residue S (sequence: GTDVKGTIFFSQEGDGPTTVT) with score 0.692.

Protein-protein interaction (PPI) analysis

It has been observed that the interacting partner proteins with query protein are IDP712 (uncharacterized LOC100283786 in Zea mays), sod3, pco095461 (Superoxide dismutase in Zea mays), GRMZM2G139680_P01 (2-cys peroxiredoxin BAS1 in Zea mays), SODA.3, 541646, GRMZM2G157018_P01 (ATP synthase subunit d, mitochondrial in Zea mays), GRMZM2G125151_P01 (uncharacterized LOC100283203 in Zea mays) and PER1 (period circadian regulator 1 in humans) with similarity scores 0.979, 0.883, 0.852, 0.847, 0.847, 0.837, 0.837, 0.835, 0.789 and 0.761 respectively (Fig. 5, Table 3). The process and function of SOD2 has also been compared with other proteins (Table 4). Other analysis such as protein interaction descriptions and other parameters are furnished in Supplemental Files (Tables S1 and S2).

Domain composition prediction

Domains play a significant role in the functional activities of proteins in the cell, and are either active regularly or during the process of evolution. The domain present in SOD2 protein sequence is Sod_Cu, which is metalloprotein that helps in the prevention of damage by free radicals with the catalysation of SOD into O² and H₂O_2. Two outliers homologous and homologous of known structures are also found (Table S3). The fragment GFHVHALGDTT was identified as the plant dismutase signature motif (42–52 residues). This was a consensus pattern in both target sequences of Zea mays (HKI 335) and template Potentilla atrosanguinea. GNAGGRUACGII signature motif was identified from heme ligand site which was conserved in both SOD2 of target and template sequences.