The amino acid sequences of Arabidopsis NPFs were used as references. The hidden Markov model (HMM) model file (PF00854) for the AtNPF gene was obtained from the Pfam database (El-Gebali et al., 2019). Then, HMME 3.0 software (Finn et al., 2015) was used to search for homologous genes in three Gossypium species (Zhu et al., 2017), with an E-value < 1e^-5, and the preliminary candidate genes were identified after we omitted incorrect and redundant members. Finally, SMART (https://smart.embl.de/smart/set_mode.cgi?NORMAL=1), PfamScan (https://www.ebi.ac.uk/Tools/pfa/pfamscan/) and the NCBI Conserved Domain Database (CDD) (https://www.ncbi.nlm.nih.gov/cdd) online websites were used to further confirm whether these candidate NPF proteins contained conserved domains. The physicochemical properties of the GhNPF proteins were predicted using the online website ExPASy (https://web.expasy.org/compute_pi/). The subcellular localizations of GhNPF proteins were predicted using the WoLF PSORT online website (https://wolfpsort.hgc.jp/). Genomic datum for Arabidopsis, G. hirsutum, G. raimondii and G. arboreum were obtained from The Arabidopsis Information Resource (TAIR) (https://www.arabidopsis.org/index.jsp) and Cotton Functional Genomics Database (CottonFGD) (https://cottonfgd.net/about/download.html), respectively. TBtools 1.098745 (Chen et al., 2020) software was used to map the locations of the genes on the chromosomes, and genes were named according to the chromosomal locations of the NPF gene family in G. hirsutum species.

To identify tandem and segmental duplication events of NPF genes, a multiple sequence alignment of full-length NPF proteins was performed by MCScanX; for this, whole-genome sequences of G. hirsutum, G. arboreum and G. raimondii and gene annotations were used. Plots of the data were created using TBtools (Chen et al., 2020) software. To assess the evolutionary constrictions on each gene pair, the non-synonymous (Ka) and synonymous (Ks) substitutions were calculated using the Simple Ka/Ks Calculator (NG) in TBtools (Chen et al., 2020). To further observe the interspecific and intraspecific homology of the NPF genes, phylogenetic trees were constructed based on the NPF protein sequences of Arabidopsis, G. hirsutum, G. arboreum and G. raimondii. The ClustalW tool of MEGA-X software (Kumar et al., 2018) was used to align the protein sequences of the cotton and Arabidopsis NPF gene family members, and then the neighbor-joining (NJ) method was used to construct a phylogenetic tree; the Poisson model was used, and the bootstrap value was 1,000. Finally, the online tool iTOL (https://itol.embl.de/upload.cgi) was used to produce a high-quality phylogenetic tree map.

The structures of the GhNPF genes were investigated on the basis of the G. hirsutum genome annotation data via the Visualize Gene Structure tool in TBtools (Chen et al., 2020). The conserved motifs of the GhNPF genes were explored via the online website MEME (https://meme-suite.org/meme/doc/meme.html), and the maximum base numbers were set to 10, with the default parameters used. The Gene Structure View tool in TBtools (Chen et al., 2020) was used to illustrate the gene structures and construct conserved motifs maps.

To understand the possible regulatory and response mechanisms of GhNPF genes, the promoter region was selected for analysis. For this purpose, the 2,000 bp nucleotide sequence upstream of the start codon of the GhNPF family members were obtained from the CottonFGD (https://cottonfgd.net/about/download.html). The online website PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) was used to screen cis-acting elements in the promoter region. The gene function of the GhNPF family in G. hirsutum was annotated with gene ontology (GO) by using DAVID bioinformatics resources (https://david.ncifcrf.gov/). ChiPlot (https://www.chiplot.online/) online analytical tools was used to plot.

To verify the expression profiles of the GhNPF genes in various tissues of G. hirsutum, the RNA-seq data for 9 tissue-specific samples of upland cotton (TM-1) (root, stem, leaf, sepal, petal, anther, pistil, ovule and fiber) and samples under salt, drought, cold and heat stress were downloaded from Zhejiang University (ZJU) (http://cotton.zju.edu.cn/) (Zhang et al., 2015). The Gossypium Resource and Network Database (GRAND) website (http://grand.cricaas.com.cn/home) was used to obtain the RNA-seq data for 9 different tissues of upland cotton (TM-1) and samples under salt, drought, cold, and heat stress from the CRI. The transcript abundance of GhNPFs in different tissues and in response to different abiotic stresses was calculated according to the fragments per kilobase of transcript per million mapped reads (FPKM) values. Heatmaps of all 98 GhNPF genes were generated using TBtools software, and Venn diagrams of candidate genes were plotted using the hiplot online website (https://hiplot-academic.com/basic/venn2).

The upland cotton cultivar Zhongmian 113 (ZM113) was grown in a greenhouse (25°C; 16 h/8 h light/darkness; humidity of approximately 60%-80%) at Gansu Agricultural University, Lanzhou, Gansu Province, China. The seeds were obtained from the CRI. Nine different organs (roots, stems, leaves, petals, sepals, anther, pistils, ovules and fibers) were collected from ZM113, which was healthy at budding and flowering stage and immediately frozen in liquid nitrogen for subsequent experiments. Healthy ZM113 plants of the same age (4 weeks old) were selected for abiotic stress treatments (heat, cold, salinity and drought). All the plants were grown in a growth chamber at 25°C before stress exposure. Each abiotic stress was applied for 0 h (control), 1 h, 3 h, 6 h, 12 h and 24 h (10 replications per treatment). Some ZM113 seedlings were subjected to cold (12°C) and heat (42°C) stress. For other ZM113 seedlings, their roots were soaked in 200 mmol/L NaCl and 15% polyethylene glycol (PEG-6000) to induce salinity and drought stresses. After the above stresses were applied, the shoot tips and young leaves were collected and immediately frozen in liquid nitrogen for subsequent experiments.

Total RNA was extracted from the shoot tips and young leaf samples collected after the stress treatments and from the tissue of nine different organs via an RNA Prep Pure Plant Kit (Tiangen, China). Two micrograms of total RNA were used to synthesize 20 µl of cDNA using FastKing gDNA Dispelling RT SuperMix (KR118) (Tiangen, China) to analyze the relative expression of the GhNPF genes in the nine organs and under the different abiotic stresses. The GhNPF gene primers used were designed using NCBI Primer-BLAST (a primer design tool) and developed by Sangon Biotech (Shanghai) Co., Ltd.; the primers used are shown in Table S1 . Real-time PCR amplification was performed using an LightCycler^® 96 Instrument together with SuperReal Premix Plus (SYBR Green) (FP209, Tiangen, China) according to the manufacturers’ instructions. The thermocycle procedure was as follows: 95°C for 3 minutes, followed by 40 cycles of 95°C for 5 seconds and 60°C for 15 seconds. All the data were normalized to those of actin (Wu et al., 2021), which served as an internal reference gene, and the relative expression of all the evaluated GhNPF genes was calculated using the 2^-ΔΔCt method (Willems et al., 2008). After normalization of the data from three independent experiments, all the data were expressed as the mean ± standard error. One-way analysis of variance (P<0.05), least significant difference (LSD) was used to evaluate the significance of each sample.

In this study, in total, 98, 52 and 51 NPF genes were verified in G. hirsutum, G. raimondii, and G. arboreum, and the G. hirsutum genes were denoted GhNPF1 to GhNPF98 according to their physical locations on the chromosome ( Figure 1 ). The details of these GhNPF gene family members and their related proteins are listed in Table S2 . The interrelated protein length (amino acids [aa]) varied greatly from 537 aa (GhNPF18) to 818 aa (GhNPF81). The predicted molecular weights (MWs) and isoelectric points (pIs) of the proteins ranged from 59,756.73 Da (GhNPF19) to 89,983 Da (GhNPF79) and from 5.38 (GhNPF52) to 9.56 (GhNPF36), respectively. With respect to the secondary structure of the GhNPF proteins, alpha-helices (Hh) and random coils (Cc) accounted for a large proportion, while extended strands (Ee) and beta turns (Tt) constituted a comparatively low proportion. Subcellular localization predictions showed that the great majority of the proteins encoded by the GhNPF genes were located at the plasma membrane, except in the cases of those encoded by GhNPF13 and GhNPF61.

The 96 GhNPF members were disproportionately located across the 26 chromosomes of G. hirsutum, and two genes (GhNPF97 and GhNPF98) were on scaffolds ( Figure 1 ). Chromosomes A03, A05 and D05 contained the greatest numbers of GhNPFs (7 members), while chromosomes A07, A13, D02, D07, D08 and D13 contained 5 GhNPFs, and they accounted for a large portion of the GhNPFs across the 26 chromosomes. In contrast, chromosomes A01, A11 and D11 contained the fewest GhNPF genes (1 member each).

To reveal the homologous locus relationships of the GhNPF gene family members between the A_t and D_t subgenomes in G. hirsutum, gene duplication events were studied using the MCScan tool to elucidate their amplification patterns. Two pairs of genes with tandem repeats were identified on chromosomes A03 and D05 (GhNPF5/6 and GhNPF63/64), respectively. In addition, 84 segmentally duplicated genes were discovered in the GhNPF gene family of G. hirsutum ( Figure 2 , Table S3 ). These results showed that segmental duplication accounted for a large proportion of the evolution of the GhNPF gene family, which reflected the dominant role of segmental repeats relative to tandem repeats in the GhNPF gene family evolution. Moreover, the intergenomic synteny analysis results between G. hirsutum and two other Gossypium species were compared to further the understand homologous gene functions and phylogenetic relationships of NPF genes ( Figure S1 ). The analysis of collinearity among the different species showed that 79 pairs of genes were collinear between G. hirsutum and G. arboreum and between G. hirsutum and G. raimondii. In conclusion, the present results provide evidence that NPF genes might undergo some genomic rearrangements during polyploidy. To better comprehend the evolutionary constraints controlling the functional divergence of the GhNPF gene family, the non-synonymous substitutions (Ka), synonymous substitutions (Ks), and non-synonymous to synonymous substitution (Ka/Ks) ratio were calculated ( Table S4 ). All duplicated GhNPF gene pairs presented a Ka/Ks ratio of <1, suggesting that the GhNPF family genes might have experienced selective pressure throughout their evolution.

To analyze the phylogenetic relationships of the NPFs among G. hirsutum, G. raimondii, G. arboreum and Arabidopsis, a phylogenetic tree comprising the NPF proteins of G. hirsutum (n=98), G. raimondii (n=52), G. arboreum (n=51) and Arabidopsis (n=53) was constructed ( Figure 3 ). The 98 GhNPF proteins clustered into three primary groups (Group I, Group II and Group III) according to bootstrap values (=1,000). There were only two GhNPF genes (GhNPF34 and GhNPF82) in G. hirsutum belonging to Group I. There were eight GhNPF genes in G. hirsutum belonging to Group II. At the same time, Group III was unevenly divided into four subgroups: III-1, III-2, III-3 and III-4. Furthermore, the GhNPF members essentially clustered into subgroups III-2, III-3 and III-4, and the number of GhNPFs in G. hirsutum was two to three times greater than that in Arabidopsis among these subgroups. Among these species, 18 pairs of paralogous genes were found—15 pairs of genes in Arabidopsis, two pairs in G. hirsutum and one pair in G. raimondii. Furthermore, 86 pairs of orthologs from G. hirsutum, G. arboreum and G. raimondii were identified, revealing the paralogous and orthologous connections among these plant species.

To research the gene structure in the evolution of the G. hirsutum gene family, the structures of the GhNPF genes were obtained by analyzing the exon/intron boundaries ( Figure 4A ). The analysis of exon/intron structure revealed relatively high structural divergence among the GhNPF genes. The number of exons in the 98 GhNPF genes ranged from three to seven, and GhNPF81 contained the most exons (n = 7). Most of the genes in Group I contained three introns, whereas in Group II, the genes contained two and four introns. Most of the genes in Group III contained three and four introns and the genes (GhNPF81) with the most introns were also included in the Group III. A total of 54.08% of all GhNPFs (53 genes) contained three introns each, suggesting that introns were gained and lost as the GhNPF gene family evolved, which might have resulted in functional diversity among the GhNPF genes.

Ten conserved motifs were detected in most GhNPF protein sequences by the use of the online MEME program, which further the similarities and differences in motif composition ( Figure 4B ). The amino acid numbers of the motifs ranged from 21 to 41 ( Figure 4C ). The number of motifs for each GhNPF was nine to fourteen ( Figures 4B, C ). Motifs 1, 2, 3, 4, 5, 6, 7 and 8 were present in all the GhNPF proteins, while motif 9 was not present only in the GhNPF18 protein of Group III. Similarly, motif 10 was not present in the GhNPF14, GhNPF62, GhNPF63, GhNPF25 and GhNPF73 proteins of Group III. In contrast, motifs 1-10 were all present in all the GhNPF members of Groups I and II. In general, almost all the GhNPF proteins within the same subgroup presented very similar motif compositions, suggesting that these GhNPF proteins have similar functions.

The cis-regulatory elements in the 2,000 bp upstream region of the 5’ end of the 98 GhNPF genes were identified and analyzed to reveal their potential response mechanisms ( Figure 5 ). We identified 55 cis-regulatory elements involved in stress responsiveness, tissue-specific expression, phytohormone responsiveness and light responsiveness. Five stress-related elements were identified, namely, DREs, LTRs, MBSs, TC-rich repeats and WUN motif–containing elements. These cis-acting elements were involved in responses to low temperature, salt stress, defense and drought. There were seven cis-acting elements associated with tissue-specific expression, namely, AREs, CAT-boxes, GC-motif, GCN4_motif, HD-Zip 1s, O2-site–and RY elements. Moreover, among these elements AREs were the most common in the GhNPF gene promoters. In addition, eleven hormone-related elements, namely, ABREs, AuxRR-core–containing elements, CGTCA motif–containing elements, GARE motif–containing elements, P-boxes, SAREs, TATC-boxes, TCA elements, TGA-boxes, TGACG motif–containing elements and TGA elements, were also found. This category included abscisic acid-responsive elements (ABREs), auxin-responsive elements (AuxRR-core–containing elements, TGA-boxes and TGA motif–containing elements), methyl jasmonate (MeJA)-responsive elements (CGTCA motif–containing elements and TGACG motif–containing elements), gibberellin-responsive elements (GARE motif–containing elements, P-boxes and TATC motif–containing elements) and salicylic acid-response elements (SAREs and TCA elements). There were also 32 cis-acting elements related to the light response, including Box 4 elements, C-boxes, G-boxes, etc. Box 4 elements and G-boxes were present in relatively high numbers within the light-responsive cis-acting regulatory elements. Interestingly, we found that the GhNPF26 gene does not contain any type of cis-acting element. Taken together, these results showed that GhNPF genes might play an important role in abiotic stress responses, defense-related signal transduction, and phytohormone responses. In addition, the genes might be involved in various light responses during G. hirsutum growth.

To analyze the expression patterns of GhNPF genes during G. hirsutum development, RNA-seq data of various G. hirsutum tissues were used in this study. The expression characteristics of all 98 GhNPF genes were determined at varying levels across different tissues and developmental stages ( Figure 6 ). The RNA-seq data of ZJU showed that 55.1% of GhNPF genes were highly expressed in vegetative organs (roots, stems and leaves) and 65.3% of GhNPF genes were highly expressed in reproductive organs (petals, sepals, anther, pistils, ovules and fibers). In addition, according to the FPKM value of CRI’s RNA-seq data, 56 out of 98 GhNPF genes were highly expressed in vegetative organs (roots, stems and leaves), and 71 out of 98 GhNPF genes were highly expressed in reproductive organs (petals, sepals, anther, pistils, ovules and fibers). A total of 52 common genes were identified in vegetative organs from two RNA-seq datasets, which verified the reliability of the data ( Figure 6C ). The expression of 14 select genes in the tissues of G. hirsutum was examined via qRT−PCR, and the results were essentially consistent with the RNA-seq data ( Figure 6D ).

To further understand the functional segregation of the identified GhNPF genes, GO was performed by DAVID based on three categories: molecular function, biological process and cellular component ( Figure 7 ). A total of 39 GhNPF genes belonged to molecular function, including transmembrane transporter activity (GO:0022857), tripeptide transporter activity (GO:0042937), dipeptide transmembrane transporter activity (GO:0071916), symporter activity (GO:0015293), low-affinity nitrate transmembrane transporter activity (GO:0050054) and nitrate transmembrane transporter activity (GO:0015112). Twenty-six GhNPF genes were involved in nitrate transport (GO:0015706), transmembrane transport (GO:0055085), nitrate assimilation (GO:0042128), oligopeptide transport (GO:0006857), dipeptide transport (GO:0042938), tripeptide transport (GO:0042939), response to nitrate (GO:0010167), response to nematode (GO:0009624), response to wounding (GO:0009611) and response to jasmonic acid (GO:0009753) in biological process. A total of 24 genes can function as an integral component of the membrane (GO:0016021) and plasma membrane (GO:0005886) in cellular component. Interestingly, some GhNPFs exist in different cell components, participate in different biological processes, and have multiple molecular functions.

To analyze the potential functions of the GhNPF genes in response to abiotic stresses, the GhNPF expression levels were evaluated via two RNA-seq datasets (those compiled by ZJU and the ICR) corresponding to plants under salt, drought, cold and heat treatments. Analysis of the expression profiles showed that six genes, namely, GhNPF5, GhNPF16, GhNPF18, GhNPF65, GhNPF69 and GhNPF86, were not expressed in any of the four treatments; moreover, GhNPF77 and GhNPF29 were not expressed in response to temperature stress, and GhNPF88 was not expressed in response to drought and salt stress ( Figure S2 ). The expression of several genes was significantly increased or decreased under the cold, heat, NaCl and PEG treatments compared with the control treatment, and the DEGs differed at different treatment time periods. According to the FPKM values of the ZJU RNA-seq dataset, 52, 36, 41 and 48 genes were highly expressed during at least two of the five different time periods (1 h, 3 h, 6 h, 12 h and 24 h) and stress treatment groups (heat, cold, salt and drought) ( Figures S3A, C ). According to the FPKM values of the CRI RNA-seq dataset, 51, 36, 51 and 31 genes were highly expressed in at least two of the five different time stress treatment groups (heat, cold, salt and drought, respectively) ( Figures S3A, C ). These results showed that after heat, cold, salt and drought treatment, 59 (heat and cold) and 38 (salt and drought) genes were upregulated, and according to the two datasets ( Figures S3B, E ), 31 genes common to the temperature (heat and cold) and saline-alkali (salt and drought) treatments could be candidate genes with resistance-related characteristics in cotton ( Figure S3F ). Among these candidate genes, the GhNPF6, GhNPF37 and GhNPF54 genes were all upregulated after heat and cold treatment, but the expression of the GhNPF64 and GhNPF27 genes was inhibited by the temperature treatments. In addition, the GhNPF28, GhNPF55 and GhNPF78 genes were all upregulated after NaCl and PEG treatment, while the GhNPF74 gene was upregulated only after NaCl treatment; PEG treatment inhibited the expression of the GhNPF39 gene. To further investigate the possible response of GhNPFs to abiotic stress conditions, by performing qRT−PCR, we analyzed the expression of 13 select genes from different tissues of G. hirsutum under different stresses ( Figure 8 ). The heat treatments (at all time points) induced the expression of GhNPF13, GhNPF31, GhNPF43, GhNPF79 and GhNPF96, and cold stress exposure at all time intervals induced the expression of GhNPF15, GhNPF43 and GhNPF56. Similarly, the salinity stress treatments induced the expression of GhNPF6, GhNPF31, GhNPF54, GhNPF57 and GhNPF96 at all time intervals, and the GhNPF31, GhNPG37, GhNPF43, GhNPF79, GhNPF96 and GhNPF98 genes exhibited increased expression in response to drought stress. Taken together, the results of our abiotic stress response gene expression analysis showed that the GhNPF gene family members in upland cotton have potential regulatory roles in the response to abiotic stress.