Genome-wide identification of 2-oxoglutarate and Fe (II)-dependent dioxygenase family genes and their expression profiling under drought and salt stress in potato

Retrieval and identification of potential 2ODDs genes from potato

Firstly, the proteome of S. tuberosum was retrieved from Spud DB (http://spuddb.uga.edu) DM1-3 v6.1 (Pham et al., 2017). Reference sequences of 2ODDs of S. lycopersicum, S. chacoense, S. melongena, Manihot esculenta, Capsicum annuum, Arabidopsis thaliana, and Nicotiana tabacum were downloaded from the National Centre for Biotechnology Institute (NCBI) protein database, which are shown in Table S1. The reference genes were Blastp (e values of <0.001) searched against the potato proteome (Verma & Singh, 2021). Secondly, Hidden Markov Model (HMM) profiles of St2ODDs were fetched from the pfam database, and the retrieved sequences were searched (e value cut-off 1e-05) against HMM profiles. The candidate sequences obtained from both methods were considered putative 2ODDs genes and were further analysed for 2ODDs domains (2OG-FeII_Oxy, pfam03171 and DIOX_N, pfam14226) using various servers like SMART (http://smart.embl-heidelberg.de/) (Letunic, Khedkar & Bork, 2021), and NCBI Conserved Domain Search (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) (Marchler-Bauer et al., 2015).

Gene structure, motif analysis and chromosomal mapping

The diversity of St2ODD genes was further analysed by studying the gene structure and conserved protein domains. To visualise the intron/exon structure, the online tool Gene Structure Display Server (GSDS) 2.0 (http://gsds.gao-lab.org/index.php) (Hu et al., 2015) was used. Conserved motifs were identified using the MEME online website (https://meme-suite.org/meme) with a maximum number of motifs, 10; optimum width of each motif, between 50 and 100 residues; other parameters were set to default values (Bailey et al., 2015). Chromosomal locations of 2ODDs were collected from Spud DB (http://spuddb.uga.edu) DM1-3 v6.1 (Pham et al., 2017), and were distributed across the 12 chromosomes of potato using TBtools (https://github.com/CJ-Chen/TBtools/releases) (Chen et al., 2020). Nomenclature was given based on the order of the chromosomal location of the genes.

Phylogenetic analysis and gene duplication

Multiple protein sequences of 205 identified genes were aligned using ClustalW in Molecular Evolutionary Genetics Analysis MEGA 7.0 (Kumar et al., 2016) (https://www.megasoftware.net/) and the phylogenetic tree was constructed with the neighbour-joining (NJ) method with 1,000 bootstrap replicates and complete deletion. The phylogenetic tree constructed was further analysed with the ITOL tool (https://itol.embl.de/upload.cgi) (Letunic & Bork, 2007). Gene duplication events were considered to have 80% or more identity with an e value 1e−10. The synonymous substitution (Ks) rates and nonsynonymous substitution (Ka) rates of duplicated 2ODD genes were calculated using Pal2nal http://www.bork.embl.de/pal2nal (Suyama, Torrents & Bork, 2006) and selection pressure was evaluated by calculating the Ka/Ks ratio (Shumayla et al., 2019).

Protein structure analysis and 3D modelling

The physio-chemical properties and subcellular localization of 2ODDs proteins were calculated using the ProtParam ExPasy server (https://web.expasy.org/protparam/) (Walker et al., 2005) and the ProtComp version 9.0 server (http://www.softberry.com) (Emanuelsson et al., 2000). SWISS-Model was used for 3D structure prediction of St2ODDs (Waterhouse et al., 2018) and visual representation of 3D structures of St2ODDs was done using UCSF Chimera (Pettersen et al., 2004). Ligand 3D models with PubChem CID 51 (2-Oxoglutaric acid), specific to FNS and GaOX, and PubChem CID 769 (Bicarbonate), specific to ACOs, were retrieved from the PubChem database in SDF format (https://pubchem.ncbi.nlm.nih.gov/) (Kim et al., 2016) and converted to PDB format using PyMOL (Schrödinger, LLC, New York, NY, USA).

Protein-ligand docking evaluation

The molecular docking of St2ODDs and their various ligands was performed using AutoDock 4. The 3D conformations of various ligands were retrieved from PubChem in SDF format: PubChem CID 51 (2-oxoglutaric acid) and PubChem CID 769 (bicarbonate). PyMOL Version 2.0 (Schrödinger, LLC, New York, NY, USA) was utilised for converting it to PDB format. The PDB files were then converted to PDBQT format using AutoDock 4 (Morris et al., 2009). Various parameters were assigned, including the addition of non-polar hydrogens and gasteiger charges. Grid boxes were generated with different dimensions in X,Y, and Z directions, as shown in Table S3. The proteins and ligands were docked with an energy range of 4 and exhaustiveness set to 8, and the best conformation was selected as having the lowest free energy (Anand et al., 2022). The protein-ligand hydrophobic and hydrogen bond interactions were represented with the BIOVIA Discovery Studio Visualizer (https://www.3ds.com/products-services/biovia/products/molecular-modeling-simulation/biovia-discovery-studio/).

Identification of abiotic stress responsive 2ODDs genes

Sequence Read Archive (SRA) data corresponding to project numbers SRP056128, SRP229183, and SRP237987 for drought, heat, and salt were downloaded from the NCBI SRA. Differentially expressed genes were identified using the Trinity-V 2.03 package. Transcript abundance was calculated in FPKM (fragments per kilobase of transcript per million mapped reads), and the heat map of expressed 2ODD genes under various stress conditions was visualised using TBtools (https://github.com/CJ-Chen/TBtools/releases) (Chen et al., 2020). Putative 2ODD genes responsive to stress were validated using RT-qPCR (Verma, Upadhyay & Singh, 2021).

Plant materials and treatments

Solanum tuberosum cv Kufri jyoti plantlets were received from the ICAR-Central Potato Research Institute (CPRI, Shimla) and maintained in a growth chamber in soil under a 16-hour (h) light/8-hour (h) dark photoperiod, 240 C, and 60% humidity. The seedlings were watered regularly. For drought stress, watering was withheld to mimic severe drought conditions for 3 days, followed by rewatering for 3 days. The leaves of drought-treated plants, after rewatering, and control plants were harvested immediately in liquid nitrogen and then stored at −80 °C. For salt stress, 4-week-old plants grown in soil were subjected to 500 mmol/L NaCl. The salt-treated plants were collected after 0, 24, 48, 72, and 96 h time courses. Collected leaf samples were stored immediately in liquid nitrogen, followed by −80 °C for further experimentation.

RNA isolation, cDNA preparation and qPCR analysis

Total RNA from the drought and salt stress samples, along with control leaf samples, was extracted from the samples (Ghawana et al., 2011). cDNA was synthesised using the Superscript III first-strand cDNA synthesis kit (Invitrogen USA) according to the manufacturer’s instructions. Elongation factor 1 alfa (ef1α) was the reference gene for expression normalization. The expression levels of St2ODD genes in potatoes subjected to drought and salt stress were measured using qRT-PCR. Primers were designed using Primer 3 software (http://primer3.ut.ee) (Untergasser et al., 2012) and the primer sequences are shown in Table S4. qRT-PCR was performed on stressed and control samples using the Bio-Rad CFX96 real-time PCR detection system. All the samples were taken in triplicate, and three technical replicates of each biological replicate were taken. PCR conditions were 95 °C for 5 min, followed by 38 cycles of 95 °C for 20 s, Tm for 20 s, and 72 °C for 20 s. Tm optimization of the stress-responsive St2ODDs was done using semi-quantitative PCR. The correlative expression data were calculated using the 2(−∆∆CT) method (Livak & Schmittgen, 2001).

Identification of 2ODDs in potato

The 2ODDs belong to the 2OGD gene superfamily and are involved in various plant metabolic pathways, which have been explored in previous studies (You et al., 2019; Farrow et al., 2014; Hagel & Facchini, 2018). Here, to identify the 2ODDs in the potato whole proteome, BLASTP searched against the reference 2ODDs of related plants, and all the putative genes were retrieved and examined for the presence of domains associated with 2ODDs, i.e., 2OG-FeII_Oxy (PF03171) and DIOX_N (PF14226), using various bioinformatics tools like SMART and NCBI Conserved Domain Search. Eventually, a total of 205 St2ODDs were confirmed to have the 2ODD domains: 2OG-FeII_Oxy (PF03171) signature domain and DIOX_N (PF14226) domain. The physiochemical properties of these identified St2ODDs were studied, and the coding sequence (CDS) lengths, genomic sequence lengths, molecular weights (MWs), isoelectric points (PIs), and grand average of hydropathicity index (GRAVY) of these genes are shown in Table S2. The MWs of St2ODDs ranged from 22.3 kDa to 108.5 kDa. St2ODD30 had the smallest amino acid sequence of 198 aa and St2ODD54 had the largest amino acid sequence of 910 aa. The exonic numbers ranged from 2 to 14 exons, and the predicted PIs values ranged from 4.68 to 8.94, suggesting that St2ODDs consist of both basic and acidic proteins. The subcellular localization analysis suggested that St2ODDs in potatoes were expressed either in the cytoplasm or secreted extracellularly. St2ODDs consist of many plant secondary metabolites that are localised in the cytosol (Xu et al., 2008). This study predicted the same and gave insight into the basic aspects of St2ODDs.

Chromosomal distribution and gene duplication of St2ODD genes

The chromosomal location of 205 St2ODDs was identified using the Spud database, and a chromosomal distribution map of St2ODDs was constructed (Fig. 1). The chromosomal distribution showed that St2ODDs are widely distributed across the 12 chromosomes of potatoes. Chromosomes 2 and 9 have the maximum number of St2ODD genes (35/205) each. Chromosome 5 has the minimum number of St2ODDs (3/205). Our results suggest that most of the St2ODDs were distributed across the proximal ends of the chromosomes. For the expansion and generation of gene families, gene duplication plays a vital role, which needed to be studied for St2ODDs in S. tuberosum. Thus, we analysed gene duplication events for 2ODDs in S. tuberosum (Hofberger et al., 2015). Gene duplication events occur due to uneven crossing over chromosome duplication. The genes with 80% or more identity located on the same chromosome were considered duplication events. Tandem duplication is defined earlier as the duplication event within the 5 Mb region of a chromosome, while other duplication events are characterised as segmental duplication (Agarwal et al., 2016). Gene duplication events were analysed for St2ODDs genes, and the identified events were confirmed for tandem duplication and segmental duplication for their genome-wide expansion. On the basis of chromosomal location, gene clusters were observed on chromosomes 1, 2, 6, 7, 9, and 11 (Fig. 1). A total of 139 duplicated genes were found on the S. tuberosum genome, which showed that 67 percent of St2ODDs were originated by duplication events, i.e., tandem and segmental (Fig. 1). Further, fifty-six percent of duplication events were tandem duplications, and forty-four percent of duplication events were segmental duplications. This suggests that most of the expansion of St2ODDs in S. tuberosum is due to tandem duplication. The synonymous (Ks) and non-synonymous (Ka) values were calculated, and their ratios (Ka/Ks) depicted the nature of selection. The Ka/Ks ratio >1 depicts positive selection, Ka/Ks = 1 depicts neutral selection, and Ka/Ks <1 depicts purifying selection or negative selection. For St2ODDs, the Ka/Ks ratios were calculated and then evaluated to determine whether the selection pressure was the driving force for evolution or not (Liu et al., 2014). In this study, 134/139 duplicated St2ODDs experienced negative or purifying selection (Table 1), and among them, thirty-three pairs of duplication events of St2ODDs had Ka/Ks less than 0.3, suggesting their lower functional divergence during evolution. However, five pairs of duplicated St2ODDs experienced positive selection.

Phylogenetic analysis of 2ODDs in potato

To understand the phylogenetic relationship between all the identified St2ODDs, a phylogenetic tree was constructed. Based on the phylogenetic tree, St2ODDs could be divided into eight groups (1–8) based on their functional aspects (Fig. 2), which were based on gene ontology studies and homology with characterised A. thaliana proteins. Amongst them, the largest group contained 2-oxoglutarate dioxygenases (2OGs). While 38 genes had gibberellin oxidase (GAOX) functions, Flavanone synthase (FNS) function was observed in 22 genes, and 10 genes have 1-aminocyclopropane-1-carboxylic acid oxidases (ACOs) function. Other functional groups were observed, including downy mildew resistance 6 (DMR6), senescence-related genes (SRGs), dioxygenase for auxin oxidation, and jasmonate-induced oxygenase genes.

Gene structure and motif analysis of St2ODDs

St2ODDs have a diversified gene structure, as the exonic number of the identified St2ODDs varied from two to 14 in number (Fig. 3A). St2ODD54 contained the maximum number of exons (14), according to Fig. 3. Next, the protein sequences of all St2ODDs were examined for the presence of different motifs using Multiple Expectation Maximisation for Motif Elicitation (MEME). We identified ten different conserved motifs, named motif 1 to motif 10 (Fig. 4), dispersed throughout the protein sequences. Four motifs (motifs 1, 2, 3, and 4) were distributed throughout the St2ODDs, and motifs 5–10 were present specifically in different clusters. Similar clusters showed similar distributions of motifs based on their functions. The conserved motifs and their sequences are shown in Table S5.

Further identified motifs were analysed, and motif 6 contained active site amino acid residues asparagine (N), tyrosine (Y), histidine (H) and aspartic acid (D). Motif 1 contained the active site amino acid residues histidine (H), arginine (R) and serine (S) which are conserved throughout the eight groups, showing the importance of motif 6 and motif 1 in the St2ODD enzyme family for functioning. The amino acids N, Y, R, and S are responsible for interacting with 2 OG, and the amino acids H, H, and D are responsible for interacting with Fe (II), which is a cofactor (Takehara et al., 2020). The conserved active amino acid residues suggest their potential role in interactions with their substrate and the functioning of these identified St2ODDs.

Protein-ligand interactions of St2ODDs

The 3D structures of representative members of FNS, GaOx, and ACOs in each group based on their function were modelled using SWISS-MODEL on the basis of homology. All the selected groups have vital roles in flavonoid biosynthesis, gibberellin biosynthesis, and ethylene biosynthesis, which needed to be studied for further understanding the substrate binding and functioning of these proteins, as shown in Fig. 5. The active amino acid residues for 2OG binding and iron co-factor binding with FNS were previously described (Sun et al., 2015). Similarly, the active amino acid residues responsible for 2OG and iron co-factor binding with GaOx were previously described (Takehara et al., 2020). For ACOs, the active residues were previously described and used as references for finding the active amino acid residues in St2ODDs (Zhang et al., 2004). The multiple sequence alignment of the three groups having FNS, GaOx, and ACO function was made and visualised, in which the active amino acids N, Y, R, and S responsible for interacting with 2 OG were conserved for binding with FNS and GaOx. For ACO, the binding sites for bicarbonate were also conserved, i.e., R and R. However, the co-factor binding, i.e., Fe-II, and the amino acids (H, H, and D) were conserved in the three groups, as shown in Fig. 6. These conserved residues suggest the functional importance of these residues for possible interactions with substrates and co-factors. In this study, a few St2ODDs have altered active amino acid residues, inferring their possible functional divergence. To further validate the involvement of these active amino acid residues, docking was performed. Hydrogen and hydrophobic interactions were studied between St2ODD29, St2ODD124 of GaOXs, and 2-Oxoglutaric acid (CID 51). St2ODD85, St2ODD87 of FNS, and 2-Oxoglutaric acid (CID 51). St2ODD118, St2ODD120 of ACOs and bicarbonate (CID 769) (Figs. 7A–7B) The docking results of St2ODD29 and 2-Oxoglutaric acid showed hydrogen bond interactions with R292 and hydrophobic interactions with S294, Q234, and E94. The docking results of St2ODD124 and 2-Oxoglutaric acid showed hydrogen bond interactions with R263, S265 and hydrophobic interactions with L205, S207. Interactions of St2ODD118 and bicarbonate showed hydrogen bonds with R244, Q239, and hydrophobic interactions with A238, Y162. Interactions of St2ODD120 and bicarbonate showed hydrogen bonds with R235 and hydrophobic interactions with R175, H234, and H177 (Figs. 7C–7D). The docking results of St2ODD85 and 2-Oxoglutaric acid showed hydrogen bond interactions with S317 and hydrophobic interactions with R315, Y232. The docking results of St2ODD87 and 2-Oxoglutaric acid showed hydrogen bond interactions with S317, L264, and I316 and hydrophobic interactions with R315 (Figs. 7E–7F). These results validated the role of the conserved amino acids in stabilising and interacting with the substrates responsible for leading to the plant pathways.

Expression pattern of St2ODDs under drought, salt, and heat stresses

Abiotic stresses play a crucial role in the growth and yield of the potato, and to explore the potential roles of St2ODDs against these abiotic stresses, expression levels of St2ODDs under drought, salt, and heat treatments were measured using available RNA-seq data. On the basis of log 2 fragments per kilobase of transcript per million fragments (FPKM) values, expression patterns were identified (Fig. 8). Under salt stress, 71 genes were differentially expressed: 39% (15/38) of gibberellin oxidases (GAOXs), which play a part in GA biosynthesis and the catabolism pathway, are expressed in salt stress (Hu et al., 2021) and 30% (3/10) of 1-aminocyclopropane-1-carboxylic acid oxidases (ACOs), which play a part in ethylene biosynthesis, are expressed in salt stress (Chang et al., 2019) (Fig. 7). Under heat stress, seven genes were upregulated, among them St2ODD113, which belongs to GAOXs, and fifteen genes were downregulated. St2ODD125 and St2ODD119 encode GAOXs and ACOs, respectively (Fig. 7). Under drought stress, twenty-nine genes were differentially expressed. St2ODD130, St2ODD54, and St2ODD25 were upregulated. Thirteen genes were downregulated under drought stress. St2ODD73, St2ODD10, St2ODD34, St2ODD99, St2ODD127, St2ODD17, St2ODD125, and St2ODD30 belong to the GaOx group. St2ODD115 and St2ODD56 belong to the ACOs group. St2ODD195, St2ODD189, and St2ODD193 belong to the FNS group. The expression patterns of the St2ODD gene family suggest their role in various stress conditions, and some of them are downregulating or upregulating in all three stresses or in groups suggesting their role in accordance with abiotic stress.

Effect of St2ODDs expression under drought and salt stress

To evaluate the role of St2ODDs under drought stress and salt stress, the gene expression levels of nine St2ODDs were measured. The gene expression levels of St2ODD54, St2ODD25, St2ODD130, St2ODD22, and St2ODD112 were measured under drought treatment after 3 days and in rewatering after 3 days by qRT-PCR. For salt treatment, the gene expression levels of St2ODD138, St2ODD76, St2ODD91, and St2ODD34 were measured after 24, 48, 72, and 96 h. The genes were selected based on the FPKM values; the genes with higher values were selected and validated via qRT-PCR.

The expression levels of St2ODDs under drought stress were analysed, and St2ODD130, St2ODD54, and St2ODD25 showed increased relative expression under drought stress with 2.26 fold change (FC), 1.71 fold change (FC), and 2.5 fold change (FC), respectively (Fig. 9A). However, the relative expression of St2ODD22 and St2ODD112 was downregulated in the drought condition with 0.47 FC and 0.45 FC, respectively, but showed a significant change in FC after rewatering.

Under salt stress, the expression levels of the St2ODDs were measured in 1-month-old seedlings at different time points, i.e., 24, 48, 72, and 96 h. Four genes showed increased expression levels, i.e., St2ODD138, St2ODD76, St2ODD91, and St2ODD34, throughout the salt stress, which is in accordance with the expression patterns of the available RNASeq data (Fig. 9B). St2ODD138 was upregulated with FC of 1.77 after 24 h, 3.04 FC after 48 h, 4.1 FC after 72 h, and 5.71 FC after 96 h. St2ODD76 was upregulated with FC of 1.42 after 24 h, 3.44 FC after 48 h, 5.58 FC after 72 h, and 8.66 FC after 96 h. St2ODD91 was upregulated consistently throughout, with FC of 3.61 after 24 h, 5.18 FC after 48 h, 7.6 FC after 72 h, and 16.7 FC after 96 h (Fig. 8B). St2ODD34 was upregulated with FC of 1.35 after 24 h, 1.99 FC after 48 h, 3.3 FC after 72 h, and 3.71 FC after 96 h.

The elevated relative expression levels suggest the potential involvement of St2ODDs with respect to abiotic stresses, which pose huge agricultural losses globally to the potato.