The genome databases of Gossypium arboreum (accession number: PRJNA382310) (Du et al., 2018), Gossypium raimondii (accession number: PRJNA171262) (Paterson et al., 2012), Gossypium hirsutum (accession number: PRJNA433615) (Wang M. et al., 2019) and Gossypium barbadense (accession number: PRJNA433615) (Wang M. et al., 2019) were obtained from the CottonGen database¹ (Yu et al., 2014). In addition, the genome databases of Arabidopsis thaliana (accession number: PRJNA10719), Oryza sativa (accession number: PRJNA448171), Zea mays (Schnable et al., 2009), Populus trichocarpa (Tuskan et al., 2006), Sorghum bicolor (accession number: PRJNA374837), Triticum aestivum (Mayer et al., 2014) and Glycine max (accession number: ACUP00000000) were downloaded from the phytozome database². The protein sequences of AtKEA1-AtKEA6 were used to construct the hidden markov model (HMM) of the conserved domain of a specific KEA gene family. The HMMER 3.0 program and the constructed HMM model were used to search the predicted KEAs from the above plant genomes. The default parameter of the e-value threshold was set at 1e-50. Then, the NCBI Conserved Domain Database³ and SMART database⁴ were used to confirm whether the candidate protein sequences contain the special domain of the KEA family (Letunic et al., 2015). The protein sequence length, molecular weight (Mw), isoelectric point (pI), grand average of hydropathicity (GRAVY) and subcellular localization of the identified KEA members were predicted from the ExPasy website⁵ and ProtComp 9.0⁶ (Artimo et al., 2012).

All the identified KEA protein sequences from G. arboreum, G. raimondii, G. hirsutum, G. barbadense, Arabidopsis thaliana, Oryza sativa, Zea mays, Populus trichocarpa, Sorghum bicolor, Triticum aestivum, and Glycine max were aligned by ClustalX 2.0 (Larkin et al., 2007). The phylogenetic tree was constructed using the neighbor-joining (NJ) method of MEGA 7.0 with the p-distance model and 1000 bootstrap replications (Kumar et al., 2016).

The chromosome physical locations of the KEA gene family were extracted from the genome annotation file information of G. hirsutum, G. raimondii, G. barbadense, and G. arboreum, and the positions of KEA genes on chromosomes were visualized with Map Chart 2.2 software (Voorrips, 2002). The replication gene pairs of G. hirsutum, G. raimondii, and G. arboreum were detected by MCScanX software (Wang et al., 2012), and gene replication was confirmed according to the following conditions: the coverage of the alignment sequence was ≥80% of the longer gene; the similarity of the regions on the alignment was ≥80%; tightly linked genes on the same chromosome were considered tandem duplications. Circos was adopted to plot the diagram of segmental duplication events on chromosomes (Krzywinski et al., 2009). KaKs_Calculator 2.0 software was used to calculate the non-synonymous mutation rate (Ka) and synonymous mutation rate (Ks) of KEA gene replication (Wang et al., 2010).

The GhKEA proteins were used for multiple sequence alignment by ClustalX 2.0. The exon-intron structures of upland cotton KEA genes were analyzed on the Gene Structure Display Server (GSDS 2.0⁷) (Hu et al., 2015). The conserved motifs of KEA proteins were identified by the MEME program. The optimization parameters were set as follows: size distribution, zero or once per sequence; motif count, 10; pattern width, between 6 and 50 residues (Bailey et al., 2006).

DNA sequences 2000 bp upstream of the KEA start codon (ATG) were retrieved from the G. hirsutum genome database and submitted to PlantCARE⁸ for analysis of cis-acting elements (Lescot et al., 2002).

The raw RNA-sequencing data of G. hirsutum TM-1 in different tissues were obtained from previously reported transcriptome data (accession number: PRJNA248163) (Zhang et al., 2015). TBtools was used to draw a heatmap, using row-scale and zero to one scale methods, showing the expression patterns of GhKEAs (Chen et al., 2020).

In this study, G. hirsutum Texas Marker-1 (TM-1) was cultivated by hydroponics with hoagland nutrient solution by Solarbio Biology Co., Ltd., and grown in a climate-controlled chamber with a light/dark cycle of 16 h at 28°C/8 h at 22°C. When the third true leaf was unfolded (approximately 4 weeks), it was treated with NaCl (300 mMol/L), KCl (0.03 mMol/L), PEG6000 (30%) and control group. Samples were taken at 0, 1, 3, 6, 12, and 24 h respectively. All the samples were immediately frozen in liquid nitrogen and stored at –80°C.

The GhKEA4 (Ghir_D12G011600.1) and GhKEA12 (Ghir_A06G009360.1) fragments of 300 nt were introduced by primers, respectively. The fragments of the above genes were then ligated into the pYL156 vector. The primers used for vector construction are listed in Table 1. The recombinant vector was transformed into Agrobacterium tumefaciens LBA4404. According to the method mentioned by Gao et al. (2016), we injected LBA4404 bacterial solution carrying pYL156 (empty vector), pYL156-GhKEA4, pYL156-GhKEA12, pYL156-CLA1 (positive control) and pYL192 (helper vector) into the cotyledons of TM-1. After 24 h of dark treatment, the cotton plants were moved to a greenhouse with 12 h of light/12 h of darkness for approximately 15 days, and then treated with NaCl. Before salt treatment, 5–6 leaves of control plants and VIGS plants were taken and weighed immediately. Then the leaves were placed in a petri dish with filter paper, placed in a 28°C incubator, set to three repeats, and weighed every other hour. Water loss rate of isolated leaves = (leaf fresh weight-leaf dry weight)/leaf fresh weight × 100%. The experiment was repeated three times independently.

After 24 h of 300 mM NaCl treatment, 0.1 g of sample powder mixed by at least 20 cotton plants were taken to determine the total chlorophyll content, soluble sugar content and proline (Pro) content, as well as catalase (CAT) activity and peroxidase activity (POD). And three biological replicates were performed. Chlorophyll was extracted with a ratio of ethanol to acetone (1:1), and the absorbance was measured at 663 and 645 nm (Sun et al., 2013). Other indicators used the kit developed by Solarbio Biology Co., Ltd., and the specific operation steps were guided according to the operating instructions.

First, control plants and plants silenced for the target gene were sampled before and after treatment with 300 mM NaCl (high salt), 0.2 mM KCl (low potassium) and 10 mM KCl (high potassium). The samples were then quickly put at 105°C for 30 min to kill and dried at 80°C for 48 h until the weight was unchanged. 0.5 g of plant sample was weighed and used for determination of Na⁺ and K⁺ contents (Bao, 2005).

Total RNA was isolated from the collected samples using the RNA-prep Pure Plant Kit (TIANGEN, Beijing, China). One microgram of RNA was reverse transcribed into cDNA using the Prime Script RT Reagent kit (Takara, Japan), and the system was diluted fivefold upon completion of reverse transcription for the next experiment. SYBR Premix Ex Taq (Takara) and the ABI 7500 Real-time PCR system (Applied Biosystems) were used to carry out quantificational real-time polymerase chain reaction (qRT-PCR) experiments. The protocol was performed as follows: step 1: 95°C for 30 s; step 2: 40 cycles of 95°C for 5 s, and 60°C for 34 s; and step 3: melting curve analysis. For each sample, three biological repeats were performed to obtain reliable results (Livak and Schmittgen, 2001). The specificity of the qRT-PCR primers was demonstrated by ePCR and melting curves. The cotton histone-3 gene (GenBank accession number AF024716) was used as an internal reference gene to normalize the expression level of the target gene (Liu et al., 2017). The data were calculated according to the 2^–ΔΔCt formula (Livak and Schmittgen, 2001). Gene specific primers for qRT-PCR were designed by Oligo 7 (Table 2).

Based on the HMM model of the KEA specific protein conserved domain constructed by the ATKEA1–ATKEA6 protein sequences, a total of 8, 8, 15, and 16 KEA members were identified from G. arboreum, G. raimondii, G. hirsutum, and G. barbadense, respectively (Table 3). All these putative genes were detected to contain the typical Na⁺/H⁺ exchanger domain (pfam: PF00999) and some members contained the TrkA-N domain (pfam: PF02254) of the KEA gene family in CDD and SMART databases. The confirmed members of the KEA gene family were named GaKEA1 to GaKEA8, GrKEA1 to GrKEA8, GhKEA1 to GhKEA15 and GbKEA1 to GbKEA16 according to the size of the e-values screened. These putative GaKEAs encoded proteins ranging from 595 amino acids (aa) (GaKEA8) to 1209 aa (GaKEA2), while GrKEAs encoded ranging from 573 aa (GrKEA6) to 1209 aa (GrKEA2), GhKEAs encoded 121 aa (GhKEA15) to 1210 aa (GhKEA2) and GbKEAs encoded 441 aa (GbKEA16) to 1209 aa (GbKEA2 and GbKEA4). Most members of the KEA family have transmembrane domains, which means that these proteins may be located on the membrane. The location of these proteins on chromosomes and the predicted molecular weight (MW), isoelectric point (pI) and grand average of hydropathicity (GRAVY) are shown in Table 3.

The results of subcellular localization prediction indicated that the proteins in KEAI were found to be located in the chloroplast (Table 3). AtKEA1, -2, -3 in Arabidopsis belong to the KEAI clade were shown to be subcellular localized in chloroplasts (Aranda-Sicilia et al., 2012; Kunz et al., 2014; Sheng et al., 2014). Online tool prediction showed that the proteins in the KEAII clade were located in the plasma membrane (Table 3). However, KEA4, -5, -6 in Arabidopsis belong to the KEAII clade were shown to be subcellular localized to the Golgi, trans-Golgi reticulum, and the prevacuolar compartment/multivesicular bodies (Zhu et al., 2018; Wang Y. et al., 2019). This is inconsistent with the online website prediction results, and the subcellular localization results of KEA in upland cotton require further validation.

Using the same method, a total of 4, 7, 7, 4,13, and 12 members were identified from Oryza sativa, Zea mays, Populus trichocarpa, Sorghum bicolor, Triticum aestivum, and Glycine max, respectively (Table 4). To examine the evolutionary relationship of KEA proteins, an unrooted phylogenetic tree was constructed using 94 KEA protein sequences from 11 species (excluding TaKEA8, TaKEA12, TaKEA13, ZmKEA6, ZmKEA7, OsKEA4 short sequences, as they did not meet the requirement of 1000 bootstrap replicates) (Figure 1). The phylogenetic tree was divided into three main categories namely KEAIa, KEAIb, and KEAII, which were consistent with the members of the AtKEA gene family (Aranda-Sicilia et al., 2012; Chanroj et al., 2012). The 11 species had distributed members in all three classifications. Among these members, the KEAII branch was the largest group, containing 46 members, while branch KEAIb was the smallest, with only 16 members. However, in the two large categories of KEAI and KEAII, the distributions of members were basically uniform. For example, 7 and 8 members of upland cotton were distributed in branch KEAI and branch KEAII, respectively. By checking the sequence characteristics, we determined that the sequence of KEAIa was the longest (except GhKEA15 and GbKEA16), and the length of KEAII was the shortest.

To determine the evolutionary relationship of KEA genes in cotton, the number and location of genes on the chromosome were analyzed (Figure 2). Each chromosome contained only one or two KEA genes. All GaKEAs and GrKEAs were distributed on 6 chromosomes, while GhKEAs and GbKEAs were distributed on six chromosomes of A_t subgenomes and six chromosomes of D_t subgenomes. Comparing the location and number of chromosomes on which the GaKEAs were located with the At subgenomic chromosomes on which the GhKEAs and the GbKEAs were located, we found that the location and number of genes distributed on these chromosomes of the KEA genes were basically the same. This phenomenon suggested that the distributions of KEA genes in the cotton genome were relatively conservative.

Gene replication events are very important for the expansion of the gene family. In general, gene replication events include tandem repeats and segmental repeats (Cannon et al., 2004; Xu et al., 2012). In this study, tandem repeat genes were defined as adjacent homologous genes on a single chromosome, and there was no more than one intermediate gene. A total of 9, 12, and 10 gene duplication pairs were identified between the A_t and D_t subgenomes of G. hirsutum and their corresponding ancestral A and D diploid genomes, respectively (Table 5). The data showed that all members of the KEA gene family were amplified only by segmental replication, which meant that segmental replication played a key role in the evolution of the KEA gene family. The syntenic relationships of putative KEA genes among two diploid genomes (G. arboreum and G. raimondii) and subgenomes in cultivated allotetraploid (G. hirsutum) were shown in Figure 3. The results showed that GaKEAs and GrKEAs were distributed among 5 and 6 chromosomes of A geneme and D genome, respectively, whereas GhKEAs were distributed among 6 and 6 chromosomes of A_t and D_t subgenomes, respectively. The Ka:Ks ratio can be used to judge whether the homologous gene is under positive selection pressure (Ka:Ks > 1) or purification selection pressure (Ka:Ks < 1). The results showed that the Ka:Ks ratios of GhKEA gene pairs were less than 1, indicating that KEA homologous gene pairs had undergone purifying selection during evolution and may have similar functions (Table 5).

Phylogenetic analysis showed that the GhKEA gene family was divided into three groups, containing 8, 2, and 5 members (Figure 4A). The sequences of most members of each subfamily were similar, indicating their close evolutionary relationship. The gene length of GhKEA15 was the shortest and that of GhKEA13 was the longest. The exon numbers of GhKEAs were ranged from 4 to 20, but most of the genes contained at least 15 exons, except that GhKEA15 contained only 4 exons (Figure 4B). The GhKEA proteins were further analyzed by the MEME program and 10 conserved motifs were identified (Figure 4C). Most GhKEA members contained multiple motifs, except GhKEA15 (1 motif), and motif 1, motif 2, motif 3, motif 4, motif 7, and motif 10 were widely distributed in these members. Motif 5 existed in the KEAIb and KEAII classes, motif 6 only did not exist in KEAIb, and motif 8 and motif 9 existed only in the KEAII class. Members of the same subfamily have similar motif characteristics, exon-intron structures and gene lengths, supporting a close evolutionary relationship.

To analyze the cis-elements that may be involved in the regulation of GhKEAs, the upstream 2000 bp sequence from the start codon (ATG) of each GhKEA gene was extracted for analysis. The cis-elements were classified into hormone response elements, stress response elements and plant growth and development elements (Figure 5). The hormone response elements included mainly salicylic acid (SA), methyl jasmonate (MeJA), gibberellin (GA), auxin (IAA) and abscisic acid (ABA). Most of the GhKEAs promoter regions contained 2–4 hormone response elements, except GhKEA11 (1 ABA response element). Of these, GhKEA5, GhKEA6, GhKEA8, and GhKEA12 contained the most hormone response elements, while GhKEA14 contained the largest number. Among these hormone response elements located in the promoter regions of GhKEAs, the largest number is the MeJA response element, followed by the ABA response element (Figure 5A). Previous studies have shown that MeJA and ABA were the main plant hormone signaling molecules under stress (Mantyla et al., 1995; Pichersky and Gershenzon, 2002), so we speculated that GhKEAs might be involved in various stress responses of upland cotton.

The stress response elements and plant growth and development elements of the GhKEAs were shown in Figure 5B. GhKEA1 contained most kinds of elements, including zein metabolism regulation elements (O2-site), low temperature response elements (LTR), endosperm expression regulatory elements (GCN4_motif), drought response elements (MBS), defense and stress elements (TC-rich repeats) and anaerobic induction elements (ARE). In addition, GhKEA7 contained the largest number of components, while GhKEA9 contained only one ARE element. Among the elements, the content of anaerobic inducible elements (AREs) was the highest, followed by Zein metabolic regulatory elements (O2 sites). In previous study, zein and its hydrolyzates have antioxidant activity (Diaz-Gomez et al., 2018), so we inferred that GhKEAs could have high antioxidant activity and antioxidant stability. These results further indicated that the KEA genes might be involved in the stress response of upland cotton.

To determine the expressions of GhKEAs in various tissues of upland cotton, the transcription data of GhKEAs in different tissues (anther, filament, pistil, petal, root, leaf, and stem) in TM-1 were analyzed. As shown in Figure 6, GhKEAs were widely expressed in different tissues, and the same gene was highly expressed in several different tissues. The expressions of GhKEAs in these tissues could be divided into three groups. Group a contained GhKEA1, GhKEA3, GhKEA5, and GhKEA6, which were expressed in anthers, petals and leaves. Group b consisted of GhKEA9, GhKEA10, GhKEA11, and GhKEA13, which were highly expressed in the pistils. The last group c was composed of GhKEA2, GhKEA4, GhKEA7, GhKEA8, GhKEA12, and GhKEA14, these genes were highly expressed in the pistils, roots and stems. The multiple expression patterns indicated that the functions of GhKEAs had been differentiated in long-term evolution.

According to the analysis of cis-elements in the promoter region and previous studies on the KEA genes in other plants, GhKEAs might be involved in the stress response. To test this hypothesis, we used the available transcriptome data to analyze the expression profiles of 15 GhKEAs under salt and drought treatments, and further verified by qRT-PCR experiment. As shown in Figure 7A, GhKEAs were regulated by PEG treatments. The transcriptome data showed that the expression levels of GhKEA2, GhKEA4, GhKEA9, GhKEA10, GhKEA12, GhKEA13, and GhKEA14 were significantly upregulated under PEG stress. However, the expression levels of GhKEA1, GhKEA3, GhKEA5, GhKEA6, and GhKEA7 were significantly downregulated. In addition, the expression levels of GhKEA8 and GhKEA11 were decreased at first and then increased. Then, we selected five genes from each of the above upregulated and downregulated groups to carry out qRT-PCR experiments and further verify their response to drought stress. The data showed that GhKEAs expressions could be regulated by PEG treatment (Figure 7B). After PEG treatment, the expressions of GhKEA2, GhKEA4, GhKEA10, GhKEA12, and GhKEA14 were increased, while the expressions of GhKEA1, GhKEA3, GhKEA5, GhKEA6, and GhKEA7 were decreased, which was basically consistent with the results in the transcriptome database.

At the same time, we also analyzed the gene expression pattern under salt conditions in publicly available RNA-seq data, and the results showed that the expressions of GhKEAs were induced by salt stress (Figure 8A). The expressions of GhKEA9, GhKEA10, and GhKEA12 were upregulated, while the expressions of GhKEA5, GhKEA6, and GhKEA11 were downregulated. What’s more, the expressions of GhKEA7, GhKEA8, GhKEA13, and GhKEA14 were upregulated at first and then downregulated, and reached the highest level at 1 h after treatment, while the expressions of GhKEA1, GhKEA2, GhKEA3, and GhKEA4 were also upregulated and then downregulated. The difference was that the expressions of these 4 genes reached the highest level at 3–6 h. The inconsistent expression patterns of GhKEAs under the same stress may be due to their different promoter elements, resulting in their possible regulation by different upstream genes and thus affecting their expression patterns. We selected 10 GhKEAs from the above groups to further investigate the effect of salt stress on their expressions by qRT-PCR (Figure 8B). The results were basically consistent with the RNA-seq data; the expression levels of GhKEA2, GhKEA7, GhKEA8, GhKEA13, and GhKEA14 were upregulated at first and then downregulated, GhKEA5 and GhKEA11 were downregulated, and GhKEA9, GhKEA10, and GhKEA12 were upregulated.

To verify the potential roles of GhKEAs in potassium absorption and transport, qRT-PCR experiments were used to observe whether the GhKEAs responded to low potassium treatment. The results are shown in Figure 9. The expressions of GhKEAs responded to low potassium treatment, and 15 GhKEAs were divided into three groups according to their expression characteristics. In group a, the expressions of GhKEAs were downregulated after low potassium treatment. The expressions of GhKEA5, GhKEA6, and GhKEA15 in group b were upregulated at first, then downregulated, and then upregulated. The gene expressions in group c first increased and then decreased. Among them, the expressions of GhKEA7, GhKEA8, GhKEA9, and GhKEA12 reached the highest level at 1 h after low potassium treatment, while the expressions of GhKEA2, GhKEA4, and GhKEA11 reached the highest level at 12 h, and only GhKEA10 reached the highest level at 3 h. These results indicated that GhKEAs might regulate the abiotic stress response to low potassium ions in the environment and participate in the absorption and transport of potassium ions in cotton.

Analysis of the promoter regions of GhKEAs revealed that the promoter regions of both GhKEA4 and GhKEA12 contain ABA and MeJA hormone response elements as well as antioxidant response elements ARE. GhKEA12 was significantly differentially expressed under low potassium stress and significantly up-regulated by drought and salt stress. The expression of GhKEA4 was significantly up-regulated under drought stress, while it also responded to low potassium and high salt stress. Then, GhKEA4 and GhKEA12 were selected for the VIGS experiments and salt treatments. The albino phenotype on the pYL156-CLA1 cotton plant ensured the success of the VIGS experiment (Figure 10A). As shown in Figure 10B, the expressions of the GhKEA4 and GhKEA12 in the corresponding VIGS plants leaves were significantly lower than the expressions in the pYL156 empty vector plants, indicating that the genes had been silenced successfully. Then these plants were treated with salt to observe the phenotype. Figure 10A showed that after salt stress, cotton plants silenced for GhKEA4 and GhKEA12 showed obvious wilting compared with control cotton plants. Subsequently, the leaves of these plants were treated with drought in vitro, to calculate the water loss rate of detached leaves. The results showed that the water loss rates of cotton leaves were significantly higher than that of control plants after silencing GhKEA4 and GhKEA12 genes, indicating that the water holding capacity of leaves decreased after silencing the target gene (Figure 10E). In order to further investigate the effects of salt stress on the physiological and biochemical characteristics of plants, the chlorophyll content, proline content, soluble sugar content, peroxidase (POD) activity and catalase (CAT) activity in VIGS cotton leaves were measured after salt stress. The results showed that after salt treatment, CAT and POD in control plants were extremely significantly and significantly higher than those in plants silenced for the target gene (Figure 10F). Furthermore, the soluble sugar, proline and total chlorophyll contents in the leaves of control plants were significantly higher than those of plants silenced for the target gene (Figures 10G–I). These results indicated that the salt tolerance of cotton plants decreased after silencing GhKEA4 and silencing GhKEA12.

In addition, the results of the expression levels of related genes encoding K and Na transporters showed that GhAKT2 was down-regulated, GhHKT1, GhPOT11, GhNHX1, and GhNHX6 were up-regulated in plants silenced for the target gene, while GhNHX2 was up-regulated in plants silenced for GhKEA4 and down-regulated in plants silenced for GhKEA12 (Figures 10C,D).

In order to preliminarily characterize the potassium ion transport activity of KEA gene, the K⁺ and Na⁺ content in cotton silenced for GhKEA4 and GhKEA12 before and after high salt, high potassium, and low potassium treatments were measured. The results showed that under control conditions (untreated), after silencing GhKEA4 in cotton plants, the Na⁺ content in stems and roots was extremely significantly higher than that in empty vector control plants, the K⁺ content in the entire plant was significantly increased, but only the Na⁺/K⁺ ratio in roots was extremely significantly increased (Figures 11B,C). After silencing the GhKEA12 gene in cotton plants, the Na⁺ content in leaves was extremely significantly higher than that in empty vector control plants, the K⁺ content in the entire plant was extremely significantly increased, and only the Na⁺/K⁺ ratio in stems was extremely significantly decreased (Figure 11A). The results showed that the Na⁺/K⁺ balance could be basically maintained after silencing GhKEA4 and GhKEA12 genes.

Under salt stress, after silencing the GhKEA4 gene in cotton plants, the Na⁺ content in leaves and roots was significantly higher than that in unloaded control plants, the K⁺ content was significantly decreased in leaves and stems and significantly increased in roots, and the Na⁺/K⁺ ratio was significantly increased in leaves and significantly decreased in roots (Figure 11), indicating that silencing the GhKEA4 gene decreased the potassium ion transport activity from the lower ground to the shoot of plants under salt stress, resulting in a higher Na⁺/K⁺ ratio in leaves. After silencing GhKEA12, the Na⁺ content in the plants was extremely significantly increased, the K⁺ content was extremely significantly increased in the roots, significantly decreased in the leaves and stems, and the Na⁺/K⁺ ratio in the plants was extremely significantly increased (Figure 11), indicating that the potassium transport activity was decreased and the sodium transport activity was increased in the plants silenced GhKEA12, resulting in the plants with a higher Na⁺/K⁺ ratio.

Under high potassium treatment, the contents of Na⁺ and K⁺ in plants silenced with GhKEA4 increased significantly, but their Na/K ratio did not change significantly (Figure 11), indicating that cotton silenced with GhKEA4 could still maintain a more stable Na⁺ and K⁺ balance. Additionally, Na⁺/K⁺ ratio was significantly increased in leaves and stems and significantly decreased in roots in GhKEA12-silenced plants, indicating that potassium ion transport activity from roots to leaves was reduced in plants with GhKEA12 gene silencing.

Under low potassium stress, the Na⁺ content in leaves and roots of plants silenced for the GhKEA4 was extremely significantly decreased and significantly increased in stems. K⁺ content was significantly increased in leaves and stems and extremely significantly decreased in roots. Na⁺/K⁺ ratio was extremely significantly decreased in leaves and stems and significantly increased in roots (Figure 11), indicating that silencing of GhKEA4 gene increased the activity of plants to transport potassium to shoots under low potassium conditions; K⁺ content in plants silenced with GhKEA12 was extremely significantly increased, the Na⁺ content in leaves and roots was extremely significantly increased in stems, Na⁺/K⁺ was extremely significantly increased in leaves, and significantly decreased in stems and roots, indicating that silencing with GhKEA12 gene improved the ability of plants to absorb potassium, but accumulated more Na⁺ in leaves, resulting in higher Na⁺/K⁺ in leaves.

Genome-wide doubling events occurring in the process of plant evolution have had a lasting and far-reaching impact on plants, and some plants have even experienced whole-genome doubling events repeatedly (Tuskan et al., 2006). Approximately 130 million years ago, the common ancestor of dicotyledons experienced a genome-wide triploid event (Jaillon et al., 2007). Then, cotton independently experienced a genome-wide pentaploid event (Paterson et al., 2012; Wang et al., 2016). Allotetraploid upland cotton enlarged the number of genes after multiple replication events. In this study, 8, 8, 15, and 16 KEA genes were identified in G. raimondii, G. arboreum, G. hirsutum, and G. barbadense, respectively. The number of KEA genes in tetraploid cotton is approximately the sum of KEA genes in G. arboreum and G. raimondii. The unbalanced distribution of KEA genes on each chromosome number proved the existence of genetic variation in the process of evolution (Paterson et al., 2012; Zhu and Li, 2013; Yu et al., 2014). Based on the classification of the KEA gene family in Arabidopsis (Maser et al., 2001), all putative KEA genes in this study can be divided into three subfamilies. The lineal homologous genes of monocotyledons tend to form lineal homologous gene pairs at the end of the branches of phylogenetic trees, while the KEA genes of dicotyledons tend to be clustered together, possibly due to the different functions of KEA proteins in monocotyledons and dicotyledons (Li et al., 2016). Furthermore, the gene members of each subfamily not only have similar gene structure, sequence length and motif structure, but also have the same results of subcellular location prediction. These results suggested that the members of the KEA gene family may show relatively conservative functions in the growth of upland cotton, especially those of the same subfamily (Palusa et al., 2007). The similarities and differences in the gene structure, domain and motif of GhKEAs may be related to the long evolutionary history and gene replication of cotton (He and Zhang, 2005). In the process of gene family evolution, tandem replication and segmented replication contributed to the emergence of gene families to a certain extent. We found that GhKEAs could be amplified only by segmental replication, indicating that segmental replication played a key role in the evolution of the GhKEA gene family. Collinear analysis showed that most of the KEA homologous gene pairs between the A_t and D_t subgenomes of G. hirsutum and their corresponding A and D diploid genomes were located in the collinear region. Based on these results, we speculated that whole genome replication was the main driving force for the expansion of the KEA gene from diploid to allotetraploid.

Maintaining the homeostasis of intracellular ions is not only the basic cellular activity needed for plant growth, but also the basis for regulating plant growth and development and coping with environmental stress (Yang et al., 2019). Previous studies have shown that AtKEA1, -3, and -4 are induced by low potassium stress (Zheng et al., 2013) and that the AtKEA gene family plays a key role in K⁺ homeostasis and osmoregulation (Zheng et al., 2013; Kunz et al., 2014; Zhu et al., 2018). In this report, the expression characteristics of GhKEAs under low potassium stress were characterized for the first time, and most of the GhKEAs were induced by low potassium stress. For example, GhKEA12, as an ortholog of AtKEA4, showed significant changes in its expression levels under low potassium stress. Many studies have demonstrated that the KEA family with Na⁺/H⁺ exchanger domains and NAD-binding (KTN) domains participates in the absorption and transport of potassium ions (Qiu et al., 2003; Bölter et al., 2020). Therefore, to preliminarily explore the involvement of GhKEAs in K⁺ transport, the contents of K⁺ and Na⁺ in the leaves of silenced GhKEA4 cotton plants and silenced GhKEA12 cotton plants were examined. In addition, silencing GhKEA4 and GhKEA12 decreased the activity of transporting potassium ions in leaves under high salt condition, resulting in higher Na/K ratio in leaves. Silencing GhKEA12 inhibited K⁺ transport activity and increased Na⁺ content in cotton leaves under high potassium stress. Silencing GhKEA12 increased K⁺ transport activity and Na⁺ absorption capacity under low potassium stress. However, the effect of silencing GhKEA4 gene on K⁺ transport activity in plant leaves was very low under high and low potassium stresses. GhKEA4 is an ortholog of AtKEA2, and it has been shown that AtKEA2 is involved in K⁺ homeostasis in chloroplasts or plastids (Aranda-Sicilia et al., 2012), which is consistent with our results. Alternatively, it has also been shown that AtKEA2 has an important function in maintaining local osmotic pressure, ionic and pH homeostasis, and the formation of thylakoid membranes (Kunz et al., 2014; Ali, 2016; Aranda-Sicilia et al., 2016). While GhKEA12 is an ortholog of AtKEA6, AtKEA6 likewise plays a role in maintaining ion homeostasis (Zhu et al., 2018), as indicated by the K⁺ uptake system (Tsujii et al., 2019). Overall, our results demonstrated that the GhKEAs were involved in regulating the dynamic balance of intracellular K⁺ during the growth and development of cotton.

It has been shown that AKT1 is involved in K⁺ uptake in the micromolar concentration range (Spalding et al., 1999); HKT1 and KUP7 have an important role in regulating K homeostasis in plants (Han et al., 2016; Wang et al., 2018); and Na⁺ and H⁺ exchange rates are significantly increased in vacuoles of plants overexpressing the AtNHX1 gene (Apse et al., 1999). We found significant changes in the expression levels of these genes in cotton silenced for GhKEA4 or GhKEA12 genes. These results suggest that there may be a “genetic compensation” mechanism for potassium ion absorption and transport in cotton.

We analyzed the cis-acting elements in the promoter region of the upland cotton KEA gene family and found a variety of stress-responsive cis-acting elements such as TC-rich repeats and MBS. Studies have shown that ABA and MeJA played an important role in regulating plant stress responses (Reyes and Chua, 2007; Tavallali and Karimi, 2019). ABA response elements (ABREs) and MeJA response elements (TGACG-elements and CGTCA-motifs) were obtained in the cis-acting elements of the GhKEA gene family. Moreover, qRT-PCR results showed that the GhKEA gene family responded positively to stress in the early or late stages of salt and drought treatments. When plants are subjected to salt stress, a large number of reactive oxygen species accumulate, leading to lipid peroxidation and interfering with the normal physiological process (Jithesh et al., 2006). POD and CAT are important antioxidant enzymes for scavenging reactive oxygen species in plants (Harb et al., 2010). At the same time, proline and soluble sugar, as important substances regulating plant cell osmotic potential, promote the scavenging of intracellular reactive oxygen species to some extent (Ullah et al., 2018). In this study, GhKEA4 and GhKEA12, which were sensitive to salt treatment, were selected for VIGS experiments. The results showed that the water loss rate of detached leaves, the contents of proline and soluble sugar, and the activities of POD and CAT in VIGS plants were lower than those in the blank control to varying degrees. AtKEA2 and AtKEA6, paralogs of GhKEA4 and GhKEA12, have been reported to confer tolerance to high Na⁺ stress in Arabidopsis (Aranda-Sicilia et al., 2012; Wang Y. et al., 2019). It has been reported that salt stress-induced production of ROS can promote K⁺ entry into the cytoplasm, thereby reducing the Na⁺/K⁺ ratio (Ma et al., 2012). We speculated that excessive accumulation of ROS in plants under salt stress may further promote GhKEAs to transport K⁺, reduce the Na⁺/K⁺ ratio, and then promote the accumulation of osmoregulatory substances to improve tolerance to abiotic stress, and the conclusion needs further verification.