The G. hirsutum cultivar GX100-2 was used for qRT-PCR and salt tolerance assay, and Zhong79 was used for qRT-PCR and drought tolerance assay. The 5-day-old cotton seedlings with the same growth state were transferred into the Hoagland liquid medium (Zhang et al., 2011) which was continuously aerated at a temperature regime of 28/20°C with 16-h light/8-h dark cycle. At the trefoil stage of the seedlings, the expression profile of candidate genes was determined by treating half of the seedlings with salt (150 mM NaCl) and the remaining half with deionized water to serve as the control. The fresh leaves were collected at 0, 1, 3, 12, and 48 h after salt stress (Zhang et al., 2011; Su et al., 2020). These samples were frozen in liquid nitrogen immediately and stored at –80°C for RNA isolation. Three biological repeats in each treatment were performed.

For the VIGS experiment, cotton seeds were soaked overnight in distilled water until the radicle sprouted. Seven sprouted seeds were planted in small pots filled with 1:1 (v/v) of vermiculate and nutritional soil and kept in the greenhouse at 28°C under a 16-h/8-h light/dark photoperiod.

The Arabidopsis thaliana Columbia ecotype (Col-0) was used as the WT.

The public expression profile of leaf under salt stress across the time course (0, 1, 3, 6, and 12 h) of G. hirsutum TM-1 was obtained from the study by Zhang et al. (2015). The expression data were gene-wise normalized, and the expression patterns were illustrated using the MultiExperiment Viewer (MeV) software.

Total RNA was extracted by hexadecyl trimethyl ammonium bromide (CTAB) and precipitated by the ammonium acetate method (Zhao et al., 2012). The cDNA was synthesized using PrimeScript™ RT Reagent Synthesis Kit (TaKaRa, Dalian, China). The gene-specific primer pairs were designed by PrimerPremier software (version 5.0) based on the coding sequences of GhAQPs and stress-related genes (Supplementary Table S1). GhUBQ7 and AtUBQ were used as internal references in upland cotton and Arabidopsis, respectively (Supplementary Table S1). Gene expression was calculated with the 2^–ΔΔCt method (Livak and Schmittgen, 2001). Each sample was analyzed with three technical replicates within each of the three biological duplicates.

Cotton leaf crumple virus (CLCrV) and tobacco rattle virus (TRV) vectors were used in the VIGS experiment under salt and osmotic stress conditions, respectively. The pCLCrV-fused cDNA fragment of magnesium chelatase subunit I (GhChlI) and TRV-fused cDNA fragment of chloroplasts alterados 1 gene (GhCLA) were used as a positive control to monitor the efficiency of VIGS experiments. The fragments targeting the candidate genes containing different recognition sites were amplified as a template and integrated into pCLCrVA or TRV2.

For the construction of the 35S:GFP-GhAQP vector, the polymerase chain reaction (PCR) product was ligated into the BamHI site and XbaI of the pCAMBIA1300-eGFP vector driven by the cauliflower mosaic virus 35S promoter. This vector was transformed into the Agrobacterium tumefaciens strain GV3101.

For the OE study, the 35S:GhTIP2;1 (pCAMBIA1300-GhTIP2;1-eGFP) vector was constructed by digesting the GhTIP2;1 coding sequence with KpnI and XbaI. The digested sequence was then inserted into a pCAMBIA1300 vector fused with green fluorescent protein (GFP) tagging, which contained hygromycin- and kanamycin-resistant genes. This vector was transformed into the A. tumefaciens strain GV3101. All the primers used in the vector construction are listed in Supplementary Table S2.

The 35S:GhTIP2;1 vector was transformed into A. thaliana (Col-0) by the floral dip method (Clough and Bent, 1998). Positive transformants were selected on the MS medium with 25 mg/L hygromycin and grew until maturation.

Subcellular localization of AQPs was predicted in WoLFPSORT¹ (Horton et al., 2007). The 35S:GFP-GhAQP vector was transformed into leaves of tobacco (Nicotiana benthamiana) K329 cultivar for the subcellular localization analysis. The signal of GFP was observed under a laser confocal scanning microscope (LSM 880, Zeiss, Germany). The 35S-mCherry-OsTIP1;1 was used as a plant vacuolar maker for the colocalization experiment (Cao et al., 2020). Excitation wavelength used in 488 nm for GFP, and the wavelength range of captured light at 515–555 nm. The excitation wavelength and gain wavelength of mCherry were 561 nm and 580–630 nm, respectively.

To investigate the promoter activity in different tissues, the 8-day-old seedlings that transferred ProGhTIP2;1:GUS were used for GUS staining. For stress treatments, the 2-week-old seedlings that were transformed into ProGhTIP2;1:GUS were treated in ¹/₂ MS medium that was supplemented with or without 10% polyethylene glycol 6000 (PEG6000), 20% PEG6000, and 150 mM NaCl for 24 h. GUS Staining Kit (Biosharp company) was employed for the GUS staining. The samples were immersed in GUS histochemical staining buffers and subsequently incubated at 37°C overnight. The samples were decolorized in 75% ethanol until the color of the negative control plants turned white. GUS activity was estimated based on the presence of blue.

All constructed vectors were transformed into A. tumefaciens strain EHA105 by a heat-shock method. The EHA105 lines contained pCLCrVA, pCLCrVA-genes, and TRV, TRV-genes, vectors were mixed with an equal volume of A. tumefaciens containing pCLCrVB and TRV1, respectively, and the mixed solution was used to infiltrate plants. The quantitative real-time PCR (qRT-PCR) was performed to further confirm that candidate genes had been silenced in VIGS experiments. The primers used in the qRT-PCR analysis are listed in Supplementary Table S1, and the primers used in the VIGS experiments are listed in Supplementary Table S2.

The cotyledons of 1-week-old cotton seedlings were infiltrated with the solution containing A. tumefaciens of pCLCrV-genes or TRV-genes according to the previous description (Gu et al., 2014; Long et al., 2020). The VIGS experiments were repeated at least three times with more than three individual plants were included.

For salt tolerance assay, plants that infiltrated with a solution containing A. tumefaciens of pCLCrV-GhPIPs after 10 days were watered by 400 mM NaCl solution regularly after every 4 days until the phenotype appeared (Long et al., 2019).

For drought tolerance assay, plants silencing of GhTIP1;1 (TRV:GhTIP1;1), GhTIP2;1 (TRV:GhTIP2;1), GhPIP1;2 (TRV:GhPIP1;2), and GhPIP2;3 (TRV:GhPIP1;2) were transferred to Hoagland liquid medium containing 13% PEG6000. TRV:GFP with no silencing fragment was used as a control. The phenotype was observed after 13% PEG6000 treatment for 1 week, and the samples of leaves and roots were evaluated for the relative water content (RWC).

To identify the stress tolerance of GhTIP2;1 overexpressed Arabidopsis (OE1 and OE3) and WT, seeds were sterilized with 5% (v/v) sodium hypochlorite and cultured on MS media, vernalized for 2 days at 4°C and incubated in a growth room (22°C, 16-h light/8-h dark cycle). Seedlings that had no significant difference in the length of primary roots were transferred to MS media with or without 150 mM NaCl and 15% PEG6000, respectively, for stress analysis. The experiments were carried out with three biological replicates, and each replicate represents 20 seedlings for each line.

To identify the chlorophyll damage caused by salt stress, leaves (infiltrated with a solution containing A. tumefaciens of pCLCrV-GhPIPs after 10 days) of the same size and from the same position were picked and washed with distilled water and then were floated in salt solution (400 mM NaCl) with the abaxial surface down. Photographs were taken after the phenotype developed. Total chlorophyll content was calculated according to the formula described in the study by Arnon (1949).

The leaves of the gene-silenced plants or seedlings of Arabidopsis after NaCl or PEG treatments were used for identifying malondialdehyde (MDA), superoxide dismutase (SOD), peroxidase (POD) activities, H₂O₂, and proline content. The MDA content, H₂O₂ content, proline content, SOD, and POD activity were determined according to the study by Ullah et al. (2018). The absorbance was measured using a UV-2550 UV-vis spectrophotometer (SHIMADZU). Three biological replications were performed. The enzyme assays were performed in three biological replicates. Proline was extracted and quantified as recommended by Zhao et al. (2009).

The experiments were conducted with three biological replicates, and each replicate represents at least 12 individuals. Each graphical plot represents the results from three repeats, and the values are displayed as the mean ± SD. Statistical significance was determined using Student’s t-tests, and P-values < 0.05 were considered statistically significant.

A total of 54, 56, and 111 candidate AQP genes were predicted in G. arboreum, G. raimondii, and G. hirsutum, respectively. The gene structure, conserved motif, and duplicate genes of AQPs in A₂, D₅, and AD₁ were analyzed, which indicated that the structure and properties of AQPs were conserved in each subfamily, yet vary among subfamilies. The physical locations of the AQP genes exhibited great diversity and complexity in the genome of Gossypium. The Ka/Ks values of duplicated AQPs were less than 1, suggesting that AQP genes have undergone strong purifying selection pressure after segmental duplication and whole-genome duplication (WGD) (Supplementary File).

Most PIPs were predicted to localize on the plasma membrane. TIPs were predicted to localize on the vacuolar membrane and plasma membrane, as well as in the cytoplasm. NIPs, SIPs, and XIPs were predicted to localize on the vacuolar membrane and plasma membrane. To determine the subcellular localization of GhAQP proteins, a C-terminal GFP fusion vector containing GhAQPs driven by a 35S promoter was constructed, and the free GFP vector was used as a positive control. Transient expression in tobacco leaves was performed by agroinfiltration. The 35S:GFP fusion protein was localized in the membrane, cytoplasm, and nucleus (Figure 1). The signals of GhPIP1;2, GhPIP2;3, and GhNIP5;1 were perceived in the plasma membrane, yet the signals of GhTIP1;1 and GhTIP2;1 were captured in the vacuole membrane. The GhSIP1;3 and GhXIP1;1 were localized to the plasma membrane and cytoplasm (Figure 1). The experimental subcellular localization of AQPs was consistent with software prediction.

The expression patterns of homologous GhAQP genes in At and Dt subgenomes were similar under salt or osmotic stress (Figure 2). Most GhAQPs (genes in the red box) that were belonging to the PIP subfamily were induced rapidly and continued to be upregulated at 3 h after salt or osmotic treatment.

A few GhAQPs (genes in the blue box) were induced at 3 h after salt or osmotic stress and were continuously upregulated until 12 h. Only three gene pairs (GhPIP1;9_A/D, GhPIP2;1_A/D, and GhSIP1;2_A/D) were continuously upregulated after salt and osmotic stresses. The expression patterns of most genes under osmotic stress and salt stress were consistent, except for GhPIP1;2_A/D.

To validate the expression pattern of GhPIPs under salt and osmotic stresses, we performed qRT-PCR (Figure 3 and Supplementary Figure S2). For easy description, we used GhPIP to represent the GhPIP_A/D gene pair. Among the 21 gene pairs, most GhPIP genes were significantly induced after 150 mM NaCl treatment, while the expression of nine genes (GhPIP1;6, GhPIP1;11, GhPIP2;1, GhPIP2;11, GhPIP1;2, GhPIP1;4, GhPIP2;8, GhPIP2;10, and GhPIP2;13) was not detectable. A total of ten pairs of GhPIPs (GhPIP1;3, GhPIP1;7, GhPIP1;9, GhPIP1;10, GhPIP2;2, GhPIP2;3, GhPIP2;5, GhPIP2;6, GhPIP2;7, and GhPIP2;9) were upregulated after salt stress. Among ten upregulated GhPIP genes, five (GhPIP1;10, GhPIP2;2, GhPIP2;3, GhPIP2;5, and GhPIP2;7) showed the highest expression at 3 h, which is consistent with the RNA-seq result. Six pairs of GhPIP genes (GhPIP1;1, GhPIP 1;5, GhPIP 1;8, GhPIP 2;4, GhPIP 2;12, and GhPIP 2;14) were alternately up- and downregulated throughout time courses of the treatment; however, these genes showed sharp upregulation at 48 h. Three genes (GhPIP2;4, GhPIP2;12, and GhPIP2;14) were significantly downregulated at 12 h, while GhPIP1;8 was downregulated at 3 h. GhPIP1;1 showed minor downregulation at 1 h and stayed to 12 h but upregulated dramatically at 48 h. GhPIP1;5 genes remained downregulated from 1 to 3 h, after that, upregulated until 48 h.

We observed that four genes (GhPIP1;2, GhPIP2;3, GhTIP1;1, and GhTIP2;1) were upregulated at the early stages after drought stress. To verify their expression levels, we performed qRT-PCR (Figure 3B). The expression of these genes reached the peak value at 2 h after osmotic stress except for GhPIP1;2, suggesting that these genes functioned in osmotic stress response. All the results showed that GhPIPs genes might be involved in the abiotic stress of cotton.

Cis-regulatory sequences are linear non-coding DNA fragments that exist in front of the promoter region. Cis-regulatory elements have various functions, which depend on their types, locations, and orientations. To expound the function of GhPIPs and GhTIPs, 1,500-bp upstream sequences of GhPIPs promoter regions were extracted and used to predict cis-elements using the PlantCARE database. A total of 272 cis-regulatory elements were detected in GhPIPs and GhTIPs. Notably, 85% were core promotor elements or binding sites of DNA binding protein. In addition, 53, 23, and 10% of the motifs were involved in light response, plant hormone-responsive, and other stress-responsive, yet 14% were undefined. These results demonstrated that GhPIPs and GhTIPs have multiple roles in cotton developmental processes and abiotic stress response. In this study, we focused on five cis-regulatory elements responding to abiotic stresses, including a cis-acting regulatory element essential for the anaerobic induction (ARE), a cis-acting element involved in low-temperature responsiveness (LTR), an MYB binding site involved in drought inducibility (MBS), a motif involved in differentiation of the palisade mesophyll cells (HD-Zip 1), and a wound-responsive element (WUN-motif) (Supplementary Table S3 and Supplementary Figure S2). In general, GhPIP genes possessed at least one stress-response-related cis-element. In this study, we discovered that most GhPIP genes had an ABA-responsive element (ABRE) that participated in ABR signaling pathways under salt stress (Supplementary Figure S2).

The GhChlI-silenced plants showed a typical photobleaching phenotype in newly grown leaves, which indicated that the VIGS system was applied successfully in GX100-2 (Supplementary Figure S3A). We also examined the expression levels of three genes (GhPIP2;2, GhPIP2;3, and GhPIP2;7) in gene-silenced plants by qRT-PCR. The gene expression was significantly decreased in gene-silenced plants than that in mock (plants transformed in empty vector) (Supplementary Figure S3B).

To further evaluate the phenotype of the gene-silenced plants under salt stress, the plants were treated with 400 mM NaCl for 12 days. Under salt stress conditions, the leaves of GhPIP2;7-silenced plants wilted more seriously, and the plant height decreased than that in mock and other gene-silenced plants (GhPIP2;2 and GhPIP2;3) (Figures 4A,B). Furthermore, the leaf disks of gene-silenced plants were incubated in 0 or 400 mM NaCl solutions for 4 days (Figures 4C,D). The leaf disks of GhPIP2;2- and GhPIP2;3-silenced plants did not produce any observable differences compared to those of mock under normal conditions. However, the leaf disks of GhPIP2;7-silenced plants showed more significant browning than those from the mock (Figure 4C). The chlorophyll content of leaf disks in GhPIP2;7-silenced plants decreased significantly, but the changes of GhPIP2;2- and GhPIP2;3-silenced plants were not significant after salt stress treatment (Figure 4D). Subsequently, we tested the MDA content, POD, and SOD activities in the leaves (Figures 4E–G). Under the control condition, the SOD activity and the MDA contents were increased, while the POD activity decreased significantly in GhPIP2;7-silenced plants. Salinity greatly reduced the activity of antioxidant enzymes (SOD and POD) and increased the MDA content in GhPIP2;7-silenced plants. The silencing of GhPIP2;2 and GhPIP2;3 in cotton did not produce any observable differences in phenotype and plant height, compared to the mock under salt stress conditions (Figures 4A,B). The trends of antioxidant enzyme activity and MDA content were not consistent under salt stress conditions (Figures 4E–G). The results indicated that the silencing of GhPIP2;7 significantly decreased cotton tolerance to salt stress.

Four genes (GhTIP1;1, GhTIP2;1, GhPIP2;1, and GhPIP2;3) were selected to validate their functions in drought stress response using the VIGS experiment. We found that the silencing of GhTIP1;1 and GhTIP2;1 resulted in the wilting and yellowing of the whole plant under drought stress (Figure 5A). The gene expression was significantly decreased in the gene-silenced plants than that in mock (Figure 5B). The chlorophyll content was increased in TRV:GhTIP2;1 plants after PEG treatment, while the RWC in roots and leaves was reduced dramatically in TRV:GhTIP1;1 and TRV:GhTIP2;1 plants under drought stress (Figure 5C). The results indicated that the silencing of TIP2;1 decreased the drought tolerance in cotton.

To evaluate the function of GhTIP2;1 in response to salt and osmotic stresses, two independent homozygous lines of the T₃ generation (OE1 and OE3) were used for the subsequent physiological experiment. The OE vector pCAMBIA1300-GhTIP2;1-eGFP contains a GFP label, so the transgenic Arabidopsis lines can be identified by detecting GFP signals. The GFP signals were detected in the root of transgenic Arabidopsis, which indicated that GhTIP2;1 had been transferred into the Arabidopsis genome successfully (Figure 6A). The 1-week-old GhTIP2;1 overexpressed lines (OE1 and OE2) and WT seedlings were transferred to ¹/₂ MS medium containing 15% PEG6000 and 150 mM NaCl, and the root length was measured after 1 week. Under control conditions, there were no phenotype differences observed between WT and OE lines. Root growth was inhibited more seriously in WT than that in OE lines under PEG treatment (Figures 6B,C). The expression pattern of GhTIP2;1 was determined in transgenic plants, harboring the GhTIP2;1 promoter that could drive the expression of the GUS reporter gene. The GUS gene was strongly expressed in cotyledons, rosette leaves, and roots in the control (Figure 6D). After 10% PEG6000 treatments and 200 mM NaCl, the GUS signal was weak in the cotyledons and rosette leaves compared with the control seedlings. GUS staining was inconspicuous in the leaves when treated with 20% PEG6000. These results indicated that GhTIP2;1 was highly expressed in roots but downregulated in leaves after salt stress and osmotic stress.

In addition, we further analyzed the antioxidant enzyme activity in GhTIP2;1-OE plants and WT under salt and osmotic stresses (Figures 6E–G). To evaluate the role of GhTIP2;1 in the oxidative stress pathway, proline, MDA, and H₂O₂ contents were measured in GhTIP2;1-OE plants and WT. The results showed that less MDA and H₂O₂ accumulated in transgenic plants than that in WT under both salt and osmotic stresses. Also, the proline content was significantly higher in transgenic plants compared to WT under salt and osmotic stresses. These results showed that GhTIP2;1-OE plant enhanced osmotic tolerance compared to WT. Furthermore, the expression pattern of stress-responsive genes was determined by qRT-PCR. The transcription levels of stress-responsive genes, including AtNHX, AtLEA, and AtP5CS, showed no significant difference among OE and WT lines under normal conditions. However, the transcription of these genes in WT and OE lines was significantly induced after salt and drought treatments (Figures 6H–J). Although the transcriptional levels of stress-responsive genes were substantially higher in the OE lines than that in WT plants after NaCl treatment, the phenotype between OE and WT lines under salt stress showed no significant difference in OE1 lines. The above mentioned results indicated that overexpressed GhTIP2;1 in Arabidopsis was more tolerant to drought stress than salt stress treatment.

In this study, we analyzed the AQPs of G. hirsutum, G. raimondii, G. arboreum, and other 34 plant species and found that the number of AQPs was consistent with the total gene number of eudicots, not monocots (Supplementary Figure S1 in Supplementary File). Before the formation of angiosperms, all plant genomes experienced two whole-genome replication events, in which eudicots and monocots experienced genome tripling and replication events, respectively (Wu et al., 2020).

Upland cotton, which is a natural allopolyploid, is an excellent plant material to explore the mechanism of genome evolution and polyploidy formation. Gene duplication is an important mechanism for increasing genetic variability and creating novel genes in plants (Moore and Purugganan, 2003). Previous analyses on biotin carboxyl carrier protein (BCCP) and phospholipase C (PLC) gene evolution in Gossypium revealed that the duplicated genes evolved independently after polyploidy formation (Cui et al., 2017; Zhang et al., 2018). To further understand the duplication events, we investigated the expansion mechanism of GhAQP genes. A total of 115 duplicated gene pairs were identified, and most of those were distributed on different chromosomes (Supplementary Figures S5, S6 and Supplementary Table S3 in Supplementary File). The result demonstrated that the expansion of GhAQP genes was mainly caused by segmental duplication. The number of AQP-coding genes in G. hirsutum was approximately the sum of G. raimondii and G. arboreum, according to the WGD event in cotton evolution (Supplementary Figure S5 and Supplementary Table S3 in Supplementary File). In this study, we observed that AQP genes in the A genome and At subgenome had common ancestors, as well as in the D genome and Dt subgenome, which indicated that AQP genes were highly conserved during cotton evolution. During the long history of plant evolution, genes have been exposed to different selective pressures, including positive selection, negative selection, and purifying selection (Flagel and Wendel, 2009). The average Ka/Ks ratio of 115 GhAQP gene pairs was less than 1, which indicated that GhAQP genes experienced purifying selection during evolution (Supplementary Table S3 in Supplementary File).

The evolutionary analysis on Gossypium AQP genes showed that most of them were greatly conserved during evolution. AQPs have typical conserved NPA motifs and ar/R selectivity filter features, which are indispensable in determining the transport channel specificity (Maurel et al., 2015). All the members of the PIP, NIP, and SIP subfamily and most of the TIPs contained the same ar/R selectivity filter. PIPs showed typical NPA motifs and highly conserved ar/R selectivity filter (F-H-T-R), which are the typical water-transporting configuration (Supplementary Table S1 in Supplementary File). These two motifs were highly conserved in PIPs of A. thaliana (Johanson et al., 2001), Z. mays (Chaumont et al., 2001), S. lycopersicum (Reuscher et al., 2013), and Brassica rapa (Kayum et al., 2017). TIPs exhibited four different forms of ar/R selectivity filter, namely, GhTIP1;1 (H-I-D-V), GhTIP2;2 (H-I-S-R), GhTIP3;1 (H-I-D-R), and GhTIP5;1 (N-V-S-L), which provided evidence for the variability of TIP subfamily (Sun et al., 2016). The NIP subfamily was quite divergent in NPA motifs and ar/R selectivity filter compared to other subfamilies in upland cotton, suggesting that the substrates for transport were diverse (Perez et al., 2014). This finding suggested that the domain of the PIP subfamily was more conserved than other subfamilies. Most PIPs, SIPs, and XIPs were predicted to be positioned on the plasma membrane, suggesting that they may regulate osmotic potential and water flows across this essential plant subcellular compartment (Figure 1 and Supplementary Table S1 in Supplementary File). TIPs were mainly located on the vacuole membrane, suggesting that the TIPs may regulate cellular osmosis and water homeostasis in cotton (Figure 1 and Supplementary Table S1 in Supplementary File). The ar/R filter in the members of different subfamilies was quite divergent, indicating their divergence in solute permeability (Supplementary Table S1 in Supplementary File).

Salt and drought stresses are the major abiotic threats to plants that affect plant growth and reduce crop yield. Excess salt may become cytotoxic to the plant, leading to cell membrane destruction (Zhu, 2001). Most homologous AQP genes in At and Dt subgenomes showed the same expression pattern under salt or drought stress (Figure 2). The PIP genes played an important role in conferring abiotic stress tolerance in plants, including drought, cold, and salt (Lu et al., 2018; Wang et al., 2019). Most of the GhPIP genes or gene pairs were rapidly induced when exposed to salt stress and osmotic stress, except for GhPIP1;2_A and GhPIP1;2_A/D (Figure 2). The structures of the PIP subfamily were highly conserved, which may explain the similar biological functions in response to abiotic stress (Supplementary Figure S1 and Supplementary Table S1 in Supplementary File). It demonstrated that most PIPs are conserved in response to abiotic stress, but the functions of a few genes are differentiated.

A total of ten GhPIP genes (GhPIP1;3, GhPIP1;7, GhPIP1;9, GhPIP1;10, GhPIP2;2, GhPIP2;3, GhPIP2;5, GhPIP2;6, GhPIP2;7, and GhPIP2;9) were significantly induced (Log₂-based value > 1) after 150 mM NaCl treatment, suggesting that GhPIP genes response to salt stress extensively (Figure 3 and Supplementary Figure S1). Most GhPIP genes showed similar expression patterns between At and Dt subgenome (Supplementary Table S4; Li W. et al., 2019). In our study, the expression of most GhPIPs was increasing at early stages after salt stress, for instance, GhPIP1;10, GhPIP2;5, and GhPIP2;7 showed the high expression at 3 h after salt stress, GhPIP1;8 reached the high level at 1 h, and GhPIP2;9 reached the peak at 12 h (Figure 3 and Supplementary Figure S1). The expression of GhPIPs increased rapidly under salt stress from 12 to 24 h and reached a high expression level at 24 h except for GhPIP1;4a_At (Dt), which reached its expression peak at 48 h, which explained the mechanism of leaves wilting in a short period after salt stress (Li W. et al., 2019). The expression of GhPIP2;7 reached a high level at 3 h after salt stress and then continuously downregulated to 48 h in this study using G. hirsutum GX100-2 (leaves); the expression of GhPIP2;1 (GhPIP2;7 in our research) increased rapidly under salt stress from 12 to 24 h and reached a high level at 24 h in G. hirsutum acc. TM-1 (root) (Figure 3; Li W. et al., 2019). Our data indicated that the response speed to salt stress was strongly related to the tissues or varieties in upland cotton.

Previous reports had demonstrated that the PIP genes played a vital role in response to salt stress and could actively regulate root and leaf hydraulics in plants (Li G. et al., 2014). GhPIP2;2, GhPIP2;3, and GhPIP2;7 were significantly upregulated in response to salt stress, which was consistent with the results in Beta vulgaris (Skorupa-Kłaput et al., 2015). AtPIP2;4 and AtPIP2;5 exhibited the upregulated expression under salt stress (Feng et al., 2018). PIP2;7, which was initially referred to as Salt-Induced MIP (SIMIP), was reported to be strongly upregulated by 150 mM NaCl treatment in the 2-week-old Arabidopsis seedlings (Jang et al., 2004). The OE of PePIP2;7 enhanced salt and drought stress tolerance of Arabidopsis and yeast (Sun et al., 2021). MsPIP2;2 conferred salt tolerance by regulating antioxidant defense system-mediated reactive oxygen species (ROS) scavenging, K⁺/Na⁺ homeostasis, and stress-responsive gene expression in Arabidopsis (Li S. et al., 2019). Taken together, salt stress stimuli resulted in a wide variety of PIP gene expression patterns. GhPIP2;2 and GhPIP2;3 were the homologous genes of AtPIP2;4, which showed the same expression pattern after 150 mM NaCl treatment (Figure 3 and Supplementary Figure S2 in Supplementary File). Plants showed no obvious salt damage phenotype whether GhPIP2;2 or GhPIP2;3 were silenced (Figure 4), which may be due to the functional redundancy of homologous genes.

There was no difference in phenotype between GhPIPs-silenced plants and control under normal conditions (Figure 4). The salt injury symptoms of GhPIP2;7-silenced plants were more severe than that in mock plants, which included yellowing, slight wilting, and dwarfing (Figures 4A,B). MDA is the product of the peroxidation reaction, which indicates the degree of peroxidation of the cell membrane and the strength of the stress reaction. Antioxidant enzymes can alleviate oxidative damage caused by salt stress in plants. Under salt stress conditions, the activity of antioxidant enzymes (SOD and POD) in GhPIP2;7-silenced plants decreased dramatically, while the MDA content increased significantly, which aggravated the salt injury phenotype of plants (Figures 4E–G). We found that the HD-Zip 1 element, which was involved in the differentiation of palisade mesophyll cells, existed in GhPIP1;9_A/D and GhPIP2;7_A/D. It was reported that PIP2;1 contributed to ABA-triggering stomatal closure through open stomata (OST)1-mediated phosphorylation (Grondin et al., 2015). Combined with the cis-elements in promoter and expression profile of GhAQPs under salt stress, GhPIP2;7 may play a positive regulatory role in response to salt stress, which affects water transport by controlling mesophyll expansion. However, the roles of PIPs in response to salt tolerance in cotton still need further research.

The silencing of GhTIP1;1 and GhTIP2;1 resulted in the wilting and yellowing of the whole plant under drought stress (Figure 5A). The RWC in root and leaves was reduced in TRV:GhTIP1;1 and TRV:GhTIP2;1 plants under drought stress (Figure 5C). GhTIP2;1 may be the key gene involved in drought stress response. To further validate its roles, we overexpressed it in A. thaliana under osmotic stress. In this study, the expression pattern of GhTIP2;1 was determined by analyzing transgenic plants harboring the GhTIP2;1 promoter that could drive the expression of the GUS reporter gene. GhTIP2;1-OE individuals grow significantly better than WT in Arabidopsis under drought stress (Figure 6A).

Proline is an osmolyte that plays an important role in oxidative stress response. The accumulation of H₂O₂ in the plant cell could cause oxidative damage, while its lower concentration correlated with drought tolerance (Ullah et al., 2018). The MDA level under stress conditions was an indicator of ROS destructive effects (Sharma et al., 2012). More proline accumulation and less MDA and H₂O₂ contents in transgenic plants suggested that the OE of GhTIP2;1 reduced the sensitivity of Arabidopsis to drought stress. NHX is Na⁺/H⁺ antiporters that maintain cellular Na⁺/K⁺ and pH homeostasis (Long et al., 2020). The upregulation of ATP5CS mainly promoted the accumulation of proline (Trovato et al., 2018). To investigate the role of GhTIP2;1 in osmotic stress, we analyzed the expression of three stress-related genes (AtP5CS, AtNHX, and AtLEA). The results showed that AtP5CS, AtNHX, and AtLEA genes were upregulated in GhTIP2;1-OE Arabidopsis plants under osmotic stress. Therefore, GhTIP2;1 may enhance the osmotic tolerance by accumulating more proline and increasing the Na⁺ efflux.