For the identification of PTG families in M. indica, PTG sequences from A. thaliana, Cicer arietinum, and Cajanus cajan were used to perform a BLAST search in the NCBI (https://www.ncbi.nlm.nih.gov) and Phytozome (https://phytozome-next.jgi.doe.gov/) databases (Johnson et al., 2008; Goodstein et al., 2012). An E-value 1e⁻¹⁰ was used as the cutoff value BLAST search. After retrieving all the sequencing, we utilized SMART and pfam tools to verify the accuracy of K⁺ transporting gene family in M. indica and deleted the sequences that lacked the conserved domains. To enhance the accuracy of our analysis, we deleted genes having common start positions in the genome representing the same locus or splice variants.

The MEGA7 software with default parameters was used for multiple sequence alignment of PTG sequences of five different species, i.e., A. thaliana, M. domestica, Oryza sativa, C. arietinum, and M. indica. Moving forward after multiple sequence alignment (MSA) of 198 protein sequences of the PTG families, with the use of the IQ TREE web server (http://iqtree.cibiv.univie.ac.at/), a phylogenetic tree was built based upon the maximum likelihood method. The tree was designed and visualized by using a web-based tool ITOL (https://itol.embl.de/) (Schultz et al., 2000).

The physicochemical properties of PTG family members were predicted by using a web-based tool ProtParam (https://web.expasy.org/protparam/) with default options. Gene ID, chromosomal locations, number of nucleotides in genomic DNA and mRNA, and number of exons were identified by using the NCBI gene database (https://www.ncbi.nlm.nih.gov/gene). The subcellular localization of K⁺ transporting proteins was predicted by using a web-based tool ProtCamp 9.0 (ProtComp—Predict the sub-cellular localization for Plant proteins (softberry.com), and then, these locations were also verified by another online tool, CELLO v.2.5 (http://cello.life.nctu.edu.tw/), by using protein sequences of K⁺ transporting gene family (Garg et al., 2016; Mehanny et al., 2022).

To predict the exons and introns of the K⁺ transporting genes, genomic DNA and cDNA sequences of all the members of K⁺ transporting genes of M. indica were retrieved from the NCBI database in FASTA format. Exons and introns in these genes were predicted using an online tool named Gene Structure Display Server (GSDS) 2.0 (http://gsds.gao-lab.org/). Motifs of the K⁺ transporting gene family were identified by an online software called MEME (https://meme-suite.org/meme/tools/meme). To predict the conserved motifs, the following parameters were selected: one occurrence per sequence (oops), number of motifs of 10, motifs with at least a width of 10, and utmost width of 50. Motif logos were also made by MEME (Guo et al., 2007; Bailey et al., 2009).

To predict the conserved domains of the K⁺ transporting gene family in M. indica, the NCBI-CDD database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) was used in which all the sequences were searched against the Pfam database. These results were also verified by the SMART tool (http://smart.embl-heidelberg.de/). Then, conserved domains of the K⁺ transporting gene family were visualized by the TB tool, using the results of the Pfam database.

For determining, the cis-acting regulatory elements of the K⁺ transporting genes in M. indica, an approximately 2,000-bp promoter region was retrieved for each gene from the NCBI gene page. After the promoter sequence was retrieved, the PlantCare database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) was employed to predict the cis-acting regulatory elements in promoter regions of K⁺ transporting genes. The PlantCare results were downloaded and opened in the Excel sheets that were used to anticipate the cis-acting regulatory elements in TBtools software (Lescot et al., 2002).

The gene location (initial and terminating position) of K⁺ transporting genes present on the chromosomes was identified by the NCBI gene database. The chromosome number for each member of the PTG family was also identified from the NCBI database. The chromosomal location was visualized by a Ritchie lab tool, Phenograms (http://visualization.ritchielab.org/phenograms/plot).

For evaluation of the gene expression configuration of the K⁺ transporting genes under abiotic stress (drought stress with bio-project), genome and gff files were downloaded from the NCBI-SRA [Home—SRA—NCBI (nih.gov)] database after the genome, and gff files were retrieved from the NCBI-genome database. The galaxy server [Galaxy (usegalaxy.org)] was utilized to obtain the fragments per kilobase of transcript per million mapped reads (FPKM) values, and the expression pattern was shown as a heat map using TB tools. Drought stress was applied to 1-year-old mango plants (Chaunsa) grown in soil-filled pots (one plant/pot) under an ambient environment (12-h photoperiod, 50%–60% relative humidity, 25°C ± 1°C day/night temperature). Drought stress was imposed by watering each pot (n = 6) with 1 L of 200 g/L PEG-6000. Leaves were collected 2, 4, and 7 days after the start of treatment. To validate the expression profiles of NGS data, real-time RT-qPCR was used. For this purpose, RNA was extracted from mango leaves using the TRIzol reagent and was quantified using NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA). With the use of 1 µg of RNA and the Maxima H-minus First-Strand cDNA synthesis kit, the RNA was reverse transcribed, and cDNA was stored at −20°C. An iTaq Universal SYBR Green Super-Mix and a qRT-PCR detection equipment (CFX96 Touch RT PCR Detection System, Bio-Rad Labs, Hercules, CA, USA) were used to perform the qRT-PCR. Gene-specific primers were mapped using an online program “Oligo Calculator” at mcb.berkeley.edu/labs/krantz/tools/oligocalc.html (accessed 26 July 2022), and the specificity of primers was verified using NCBI-primer BLAST algorithm (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) (retrieved on 26 July 2022). Actin (LOC123192663) was used as a reference gene (Yao et al., 2022). The 2^−ΔΔCT method is used to calculate respective gene expression levels on the basis of three biological replicates (Livak and Schmittgen, 2001).

A total of 37 PTGs were identified in M. indica genome ( Table 1 ). Physiochemical characteristics of K⁺ transporter in M. indica, such as the total number of amino acids (aa), molecular weights (MW), aliphatic index, gravy, and hypothetical isoelectric points (pI), were identified using the web-based tool Expasy-protparam. Subcellular localization was also identified using cello 2.0. The length of protein ranged from 340aa to 1208aa. Molecular weight varied from 39.05 to 130.77 kDa, and the pI values start from 4.92 to 9.47 in which 16 of the proteins were acidic and 20 were basic. Interestingly, the MiHAK1.1 protein possessed a pI value of 7.01. GRAVY values of AKT and KAT-like proteins were negative (hydrophilic), while the rest of the proteins were hydrophobic, and subcellular localization showed that the location of these proteins is in the plasma membrane ( Table 1 ).

To predict the functional properties as well as phylogenetic relationships of the K⁺ transport-related genes, using protein sequences of PTGs of O. sativa, A. thaliana, C. arietinum, M. domestica, and M. indica, a phylogenetic tree was constructed. The members of the tree were divided into five groups: i) HAK, ii) HKT, iii) KEA, iv) TPK, and v) Shakers (AKT and KAT, and GORK and SKOR). HAK, HKT, and KEA are transporters, and TPK and Shakers are Channels. In phylogenetic analysis, some orthologous pairs (MiHAK1/MdHAK1.2 and MdHAK6.2/MiHAK6.1), co-orthologous groups (MdHAK1.1/CarHAK1/MiHAK1.1, MdHAK7.1/MdHAK7.2/MiHAK7, and MdHAK5.4/CarHAK5/MiHAK5), and paralogous groups (MiHAK8/MiHAK8.1, MiHAK2.1/MiHAK2, and MiHAK5.3/MiHAK5.1/MiHAK5.2) of M. indica were identified in the HAK family among A. thaliana, O. sativa, and M. domestica. In the HKT family, one orthologous pair was identified in A. thaliana (MiHKT1/ATHKT1). In the KEA family, two orthologous pairs (MiKEA6/CarKEA6 and MiKEA4/OsKEA4) and two co-orthologous groups (MdKEA3.1/MdKEA3.2/MiKEA3 and MdKEA5.1/CarKEA5/MiKEA5) were identified in the genome of M. indica. In the AKT and KAT sub-families, paralogous (MiAKT1/MiAKT1.1 and MiKAT2/MiKAT2.1) and orthologues groups (AtAKT6/MiAKT6, MiAKT2/MdAKT2.1, and MdKAT3.1/MdKAT3.2/MiKAT3) were identified in M. indica. In the GORK and SKOR sub-families, one orthologous pair (MiSKOR1/CarGORK) was identified; in the TPK family, some orthologous and co-orthologous groups were identified (MdTPK4.1/MdTPK4.2/MiTPK4, MiTPK5/CarTPK5, and MiTPK1.2/MdTPK1.3); some paralogous groups were also identified (MiTPK6/MiTPK6.1 and MiTPK1.1/MiTPK1) in the M. indica genome ( Figure 1 ).

To acquire an understanding of the structural properties of PTGs, we explored the intron and exon architecture of these genes in M. indica. The results of this analysis demonstrated that the number of exons varies from 2 to 20. In HKTs, there were only two exons. In HAKs, the number of exons ranged from 7 to 12. Similarly, 19–21, 2–5, and 10–13 exons were predicted in KEA, TPK, and Shaker families, respectively ( Figure 2A ). In some regions, CDS are less in concentration and illustrated in green color in the legend. Similarly, conserved motifs were predicted and visualized by using MEME suite. A total of 10 conserved motifs were predicted in all proteins. Among these, motifs 1–8 and 10 were present in all proteins. Motif 9 was detected in only KEAs, TPKs, and Shakers. This motif possesses the characteristic GYGD motif, which acts as a selectivity filter for K⁺ ( Figure 2B ).

The conserved domains present in the PTG proteins in M. indica were analyzed using NCBI-CDD. After this, further analysis was performed in phases and visualized by using the software Tbtools. All the members of PTGs possess conserved domains like (the K_trans superfamily, ions_trans_2, Na_H_Exchanger, KHA, FRQI superfamily, PLN00149, and PLN00151). Not only the most conserved domains are visualized but less conserved domains are also shown and can be visualized. The least conserved domains are from EFhs Superfamily, which revealed that some proteins function distantly ( Figure 3 ).

In M. indica, 36 of 37 genes were mapped on 14 (out of 18) chromosomes ( Supplementary Figure 1 ). The locus of MiHAK8.1 gene was not found. As of now, this gene has not been assigned to any linkage groups on scaffolds. Hence, its location is not displayed on the map. There were no PTGs on chromosomes 3, 7, 11, and 16. Chromosome 8 contains the maximum number (5) of genes, while only one gene was present on chromosomes 1, 2, 15, and 18. Almost 40% of members (13 out of 36) were present on chromosomes 8, 12, and 13 collectively ( Supplementary Figure 1 ).

Gene expression is controlled by cis-regulatory sequences, such as enhancers and promoters, which play crucial roles in modulating the development and physiology of an organism. In the current study, we have identified 60 types of cis-regulatory elements in the 2,000-bp promoter sequences of 37 PTGs in M. indica. The promoter sequence of each gene was found to be rich in these regulatory elements. Most abundant of them all were found to be TATA-Box (TATA-box is capable of defining the direction of transcription and also indicates which DNA strand should be read by the transcriptional machinery), CAAT-box (CAAT-box serves as a marker for the binding site of the transcription factor for RNA), O2-site (assists in regulating two transcription factors, those associated with the metabolism of carbon and amino acids as well as the resistance to abiotic stress), AT-rich element (a replication complex is formed at this specific site and DNA synthesis is initiated), CAT-box (cis-acting regulatory element related to meristem expression), ARE (regulatory element essential for anaerobic induction), MRE (MYB binding site involved in light responsiveness), MBS (MYB binding site involved in drought-response), BOX-4 (conserved DNA module involved in light responsiveness), TCA-element (regulatory element involved in salicylic acid responsiveness), ABRE (Abscisic Acid-Responsive Element), and P-box (gibberellin-responsive element); in addition to these, other regulatory elements were also identified, which were not present in the promoter sequence of all the genes and were found to be present in a very small number but could have a potentially significant role in the life of M. indica plant, which includes Circadian, MREG-box, LTRGCN4_motif, TCT-motif, CGTCA-motif, TGACG-motif, G-Box, GA-motif, chs-CMA2aAT1-motif, ATCT-motif, LAMP-element, GATA-motif, Box III, AE-box, LS7TCCC-motif, TC-rich repeats, chs-Unit 1, m1ATC-motif, A-box, 3-AF1 binding site, TGA-element, TATC-box, I-box, GARE-motif, GT1-motif, AuxRR-core, ACECCAAT-box, SARECAG-motif, AT-rich sequence, GC-motif, RY-element, GTGGC-motif, chs-CMA1a, Box II, Gap-box, MBSIHD-Zip 1, MSA-like, Sp1, ACA-motif, TGA-box, L-box, and WUN-motif ( Figure 4 ).

To check the gene expression profile of PTGs in M. indica leaflets and roots under drought stress, the gene expression data of M. indica were downloaded from the NCBI-SRA database ( Figure 5A ). The differentially regulated genes were either upregulated in both leaves and roots (MiKAT2.1, MiHAK2, MiTPK5, MiHKT1, MiHKT7, MiTPK1.1, MiHAK12, MiHAK1.1, MiAKT6, and MiHKT1), upregulated in leaves and downregulated in the root (MiAKT1, MiKEA2, MiAKT2, MiKEA3, MiKEA5, MiHAK8, MiHAK5.2, MiHAK5.1, MiKAT3, and MiTPK4), and upregulated in root and downregulated in leaves (MiHAK5, MiTPK6.1, MiAKT1.1, MiTPK1.2, MiHAK1, MiHAK6.1, and MiHAK2.1), while the rest of the genes were non-responsive ( Figure 5A ). Interestingly, MiTPK1.2, MiHAK1, MiHAK2.1, HAK6.1, and MiAKT1.1 were most upregulated in roots, and MiKEA2, MiAKT2, and MiAKT1 were upregulated in leaves.

RNA-seq data were also used to analyze the gene expression profile of PTGs during the five stages of fruit ripening in M. indica. SRA data for the bioProject (PRJNA797728) were retrieved from the NCBI-SRA database. Gene expression profiles were divided into two groups, i.e., those genes that were upregulated in five different stages of fruit ripening (MiAKT6, MiHAK1.1, MiKAT2, MiKAT2.1, MiHKT1, MiTPK1.1, MiHAK7, and MiHAK12) and those genes that were highly downregulated in leaves and fruit ripening according to FPKM values (MiHAK1, MiHAK5.1, MiHAK5, MiHAK5.2, and MiTPK1.2), and the rest of the genes showed non-responsive or non-significant variations of expression ( Figure 5B ).

Real-time amplification (qRT-PCR) was performed to confirm the expression of nine selected PTGs in M. indica leaves. These genes were selected based on their expression values in drought stress response ( Figure 5A ) and higher expression during fruit development ( Figure 5B ). In response to drought stress, nine genes, including MiAKT6, MiHAK1.1, MiKAT2.1, MiHKT1, MiTPK1.1, MiHAK1, MiHAK5.1, MiHAK5, and MiHAK5.2, were selected for qPCR-based quantification for 7 days, and samples were taken at the second, fourth, and seventh days. Under drought stress, MiHAK1 was regulated differently. Although according to RNA-seq analysis MiHAK1 was downregulated in response to drought stress, it was found to be upregulated up to 3.5-fold in qPCR. Similarly, MiHAK1.1, MiHAK5.2, MiTPK1.1, MiKAT2.1, and MiAKT6 were significantly upregulated after drought stress. On the contrary, MiHAK5, MiHAK5.1, and MiHKT1 were significantly downregulated ( Figure 6 ). There was a significant correlation between the results of qPCR and RNA-seq data. However, some variations were also observed. Though similar plant growth and stress conditions were applied, such variations are potentially associated with sampling time in a day, plant age, or even genotypic differences.

The current study identified 37 potential genes of the K⁺ transport system in M. indica genome. The number of PTGs is more or less similar in major plant groups represented by V. vinifera, Glycine max, Triticum aestivum, Vigna radiata, C. cajan, and A. thaliana (Cuéllar et al., 2010; Isayenkov et al., 2011; Chen et al., 2017; Azeem et al., 2018; Azeem et al., 2021a; Azeem et al., 2021b; Ge et al., 2020; Siddique et al., 2021; Cai et al., 2021). It not only indicates the evolutionary conservation of PTGs in these plants but also predicts the existence of similar genetic mechanisms for K⁺ in other fruit trees. A total of 202 PTG protein sequences from five species—A. thaliana, M. indica, M. domestica, O. sativa, and C. arietinum—were divided into six groups according to phylogenetic analyses (HAK, HKT, KEA, AKT and KAT, GORK and SKOR, and TPK). In contrast, the phylogenetic relationship analysis showed that PTGs were more closely related to MdPTGs than AtPTGs. This outcome confirmed the finding that apple and mango also showed a closer association than Arabidopsis. It is also supported by the fact that various sequence features were common in M. indica and other species. For example, the transmembrane domains of MiHKT1.1 and AtHKT1;1 are eight. Similarly, the presence of the consensus motif GVVYGDLGTSPLY is a characteristic feature of HAK transporters (Cheng et al., 2018). This motif was present with a minor modification as GVVYGDLG(I/T)SPLY. This variation of a motif is also present in HAK proteins in A. thaliana and Camellia sinensis (Yang et al., 2020). KEA members of both M. indica and A. thaliana species showed conservation of “FLLFxxGLE and GEFAFVxxxxA” motifs. Previously, it has been established that a divergence in the amino acid sequences of proteins is related to the functional divergence of proteins and vice versa (Sangar et al., 2007). Among K⁺ channels, the presence of ANK, Ion_trans_2, and KHA domains in the Shaker proteins, similar to Arabidopsis K⁺ channels, suggests their functional similarities (Keisham et al., 2018). Likewise, the occurrence of the significantly conserved residues, i.e., GYGD and RSXpS/pTXP, in MiTPKs indicates potential functional similarities in M. indica and other plants (Ge et al., 2020; Siddique et al., 2021).

As one of the most abundant cations in plant cells, K⁺ performs a dominant role in various biological processes of plants (Ragel et al., 2019). The influx/efflux of K⁺ (mediated by channel and transporter proteins) regulates its concentration in the plant body. K⁺ channels and transporters are important contributors to plant growth and development (Sharma et al., 2013). In the current study, members of HAK, HKT, TPK, and Shaker families were differentially regulated by drought stress. These results were more or less in agreement with RNA-seq and qPCR analyses. Previously, the role of potassium channels and transporters has been documented in the drought stress response of various plants (Ahmad et al., 2016a; Ahmad et al., 2016b; Chen et al., 2017; Hassan et al., 2017; Cai et al., 2019; Qi et al., 2019; Singh et al., 2021). The overexpression of the OsAKT1 and OsHAK1 significantly affects potassium nutrition and drought stress tolerance of rice (Ahmad et al., 2016a; Chen et al., 2017). Our results also proposed similar roles of K⁺ transporters. Overexpression of the potassium channel TPKb in small vacuoles confers osmotic and drought tolerance to rice (Ahmad et al., 2016b). Taken together, these findings suggest a potential role of PTGs in the stress response of the M. indica tree.

The essential nutrients (N, P, and K) promote an increase in biomass yield by improving root growth, activating cellular enzymes, ameliorating photosynthesis, conserving energy, and positively affecting other fundamental processes in the plant body. Potassium sources can influence the K⁺ release and fruit yield of M. indica (Wang et al., 2022). It emphasizes the importance of potassium uptake and transport system from root to leaves and fruit. The current study reports significant differential regulation of K⁺ transport during five stages of mango fruit development. Previously, it is reported that supplemental foliar applications of K⁺ can improve fruit quality in Cucumis melo (Lester et al., 2005), Solanum lycopersicum (Liu et al., 2021), C. sinensis (Wu et al., 2021), and M. indica (Sarker and Rahim, 2013). The findings of the current study complement the potential role of PTGs in improving mango fruit growth, development, and quality.