The complete G. hirsutum genome sequence data were obtained from COTTONGEN¹ (Yu et al., 2014). The protein sequences of five additional plant species comprising Physcomitrella patens, Selaginella moellendorffii, Oryza sativa, Theobroma cacao, and Populus trichocarpa were retrieved from the JGI Phytozome database² and Genbank database³. The amino acid sequences of KNOX proteins from Arabidopsis thaliana, which were used as query sequences to search for cotton KNOX ortholog proteins in local BLAST with BlastP (with an threshold value of E ≤ 1e-5), were accessed from TAIR 10⁴. Then, the collected KNOX-like candidate proteins were subjected to SMART for further selection based on their conserved domain⁵.

Multiple sequence alignment was performed with ClustalW⁶. The conserved KNOXI, KNOXII, and ELK domain sequences of cotton and Arabidopsis KNOX proteins were aligned. A phylogenetic tree was constructed from full-length KNOX amino acid sequences of seven plant species, using the neighbor-joining method combined with a bootstrap analysis and the Jones–Taylor–Thornton substitution model as implemented in MEGA7.0. Branch support was estimated by performing a bootstrap analysis with 1000 replicates (Tamura et al., 2007).

Chromosome size and gene location information for GhKNOX genes were extracted from the gene annotations (gff3) file accessible from the Gossypium hirsutum genome. MapChart 2.2 software was used to determine the distribution of the genes on the G. hirsutum chromosomes. The exon and intron structure was displayed using the GSDS 2.0 online server⁷. The collinearity of gene pairs in the GhKNOX family were mapped to generate a collinearity map using Circos software.

All upland cotton plants were grown in the field at the Henan Institute of Science and Technology. Different tissues were sampled from plants of the cultivar ‘TM-1.’ For stress treatments, seeds were sown in plastic pots under a 14 h light/10 h dark photoperiod at 28°C until the seedlings attained the second leaf expanded stage and were then treated with 20% polyethylene glycol 6000 (PEG) or 200 mM NaCl. Leaves were harvested at 0, 1, 3, 6, 12, and 24 h, immediately frozen in liquid nitrogen, and stored at –80°C for total RNA extraction. Shoot meristems were harvested from plants of the early maturing cultivar ‘Zao1’ and the late-maturing cultivar ‘CCRI50’ for RNA-sequencing (RNA-seq) from the fourth leaf expanded to the seventh leaf expanded stages. The raw reads were processed to retain only clean reads by removing the adaptor sequences, low-quality sequence reads (Q < 20), and poly-N stretches (>10%). The clean reads were mapped to the upland cotton reference genome to obtain unigenes using the Tophat2 software (Kim et al., 2013b). Expression of KNOX genes in different tissues and the cold, heat, salt and drought stress treatments was analyzed using raw RNA-seq data. The raw RNA-seq data were downloaded from the NCBI Sequence Read Archive⁸. The RNA-seq expression analysis was conducted using TopHat and Cufflinks. Gene expression was expressed as fragments per kilobase of transcripts per million mapped reads (FPKM). A heatmap was generated using TBtools (Chen et al., 2020).

Total RNA was isolated from samples using a plant RNA purification kit (Tiangen). The first-strand cDNA was synthesized using the PrimeScript™ 1st Strand cDNA Synthesis Kit for RT-PCR (TaKaRa). Transcript levels were determined by quantitative real-time PCR (qRT-PCR) analysis using the Q6 Real-Time PCR System (Applied Biosystems) and SYBR Premix Ex Taq (2×) (TaKaRa). To normalize these samples, GhACTIN was as an endogenous control. Determination of reaction specificities and data processing were performed as described in previous study (Schmittgen and Livak, 2008). Gene-specific primers used for the PCR are listed in Supplementary Table 1. Three biological replicates were analyzed. The significance of differences between means was determined using analysis of variance implemented in SAS software (^∗P < 0.05, ^∗∗P < 0.01). The data were graphed using GraphPad Prism 5.

The basic information for 82 early and 67 modern cultivars from a core collection of upland cotton and the relative genomic variants were downloaded from the Hebei Agricultural University website⁹. Single-nucleotide polymorphisms (SNPs) in the KNOX genes were detected, based on the genomic location of the genes, and the number of SNPs per gene was scored using Excel 2010. The fixation statistic (F_st) was calculated with Genepop 4.0 software (Rousset, 2008). The genes F_st > 0.45 were identified as putative sites under selection during improvement (Li et al., 2018).

Based on previously described VIGS assay method (Gao et al., 2011), the genes GhKNOX2-A, GhKNOX10-A, GhKNOX14-A, GhSTM2-A, and GhSTM3-A/D were amplified from the ‘CCRI50’ cDNA library and inserted in the pCLCrVA vector. Gene-specific primers used for the VIGS assay are listed in Supplementary Table 1. The recombinant vectors were transformed separately into Agrobacterium strain GV3101. The GV3101 cells harboring recombinant plasmid were mixed with cells carrying pCLCrVB (1:1 ratio). The GV3101 cells were cultured then were injected into 10-day-old cotton cotyledons. The cotton plants were analyzed with regard to their gene expression profiles and phenotypes under salt stress. The inoculated cotton plants were grown in a greenhouse at 22°C under a 16 h light/8 h dark photoperiod. The content of malondialdehyde (MDA) and activity of peroxidase (POD) were assessed using a MDA assay kit and POD Assay Kit (Nanjing Jiancheng). The analysis was repeated three times, and each data type was analyzed from a sample of at least five plants in each independent biological experiment. The significance of differences between means was determined using Student’s t-test (^∗P < 0.05, ^∗∗P < 0.01).

We identified putative KNOX genes in the reference genome of G. hirsutum. Forty-four GhKNOX genes were identified. The GhKNOX genes were named on the basis of the similarity of the encoded amino acid sequence with that of Arabidopsis orthologs; ‘A’ and ‘D’ indicated derivation in the A and D subgenomes, and ‘a’ and ‘b’ were used to distinguish the corresponding paralogs of the same Arabidopsis ortholog. Thus, the 44 putative KNOX family genes were named GhKNOX1 to GhKNOX7 and GhSTM1 to GhSTM3, and GhKNL1 was identified in a previous study (Gong et al., 2014). The other genes identified had no highly orthologous counterparts in Arabidopsis and were named GhKNOX8 to GhKNOX14. The cotton KNOX family genes encoded a peptide ranging in length from 102 to 865 amino acids, the molecular weight ranged between 11.31 and 98.44 KDa, and the isoelectric point value ranged from 4.08 to 8.81 (Table 1).

To explore the evolutionary relationships of KNOX proteins among cotton and six other plant species, a neighbor-joining tree was constructed based on a multiple alignment of KNOX amino acid sequences. The KNOX proteins were divided into KNOXI, KNOXII, and KNATM clades (Figure 1). The KNOXI clade comprised the STM, KNAT1, KNAT2, and KNAT6 homologs derived from ferns, lycophytes, and angiosperms, and 24 GhKNOX proteins were clustered in this clade. The KNOXII clade comprised KNAT3, KNAT4, KNAT5, and KNAT7 homolog proteins. GhKNOX1-A/D were clustered in the KNATM clade. Most cotton KNOX proteins showed higher similarity with proteins from cacao and poplar; these genes were consistently clustered closely on one branch in the phylogenetic tree. Based on the classification of Arabidopsis KNOX proteins, subclades I and II belonged to the class KNOXI, and clades III and IV belonged to KNOXII and KNATM, respectively (Figure 2A).

Most (30 of 44) of the G. hirsutum KNOX genes contained four introns and five exons, and eight KNOX genes contained only three introns and four exons (Figure 2B). Only GhKNOX14-D incorporated one intron and two exons, and four genes (GhKNOX1-A/D, GhKNOX5a-A, and GhKNOX6-D) included two introns and three exons. The most highly similar exon and intron structures were observed in cotton genes within the same phylogenetic clade, thus supporting the reliability of the phylogenetic analysis. A multiple alignment of protein sequences was generated to detect the KNOX domain motifs in Arabidopsis and G. hirsutum (Supplementary Figures 1, 2). Four G. hirsutum proteins (GhKNOX1-A/D, GhKNOX5a-A, and GhKNOX6-D) contained only the KNOXI and KNOXII domains and lacked the ELK domain and homeobox KN binding domain. The GhKNOX10-A protein lacked the DNA-binding domain.

Among the 44 G. hirsutum KNOX genes, 41 members were located on 20 of the 26 chromosomes assembled in the G. hirsutum genome v1.1, and the remaining three genes were located on three unmapped scaffolds (scaffold4079, scaffold503, and scaffold1256) (Figure 3). The number of KNOX genes per chromosome ranged from zero to five. Chromosomes A05 and D05 carried five genes, whereas no KNOX gene was detected on chromosomes A01/D01, A04/D04, and A09/D09. The KNOX genes located on homoeologous A and D chromosomes was conserved identical except for A02/D02, A03/D03, and A06/D06. The circos software was used to analyze GhKNOX gene duplication events in the upland cotton genome (Figure 4). The GhKNOX genes were unevenly distributed in A and D subgenomes, and specific duplications also occurred in the two subgeomes. More than ten GhKNOX genes were located in the A and D subgenome regions, respectively. Chromosomes A01/D01, A04/D04, A07, and A09/D09 did not contain any duplicated genes, whereas chromosomes A05/D05 and A06/D06 harbored the highest number of duplications. Chromosomes A03/D03 had three genes, but only one of them was paralog gene. Chromosomes D07 had one gene, while chromosomes A07 had no paralog gene. The collinearity analysis indicated that GhKNOX genes diverged from a common ancestor, but these genes were not conserved in the A and D subgenomes.

Gene expression in different tissues may be associated with diversity in biological functions. The expression patterns of GhKNOX genes in ten organs (root, stem, leaf, torus, sepal, bract, anther, filament, fiber, and ovule) were analyzed (Figure 5). Among these genes, six genes (GhKNOX1-A, GhKNL1-D, GhKNOX3a-A/D, GhKNOX3b-D, and GhKNOX5b-A) in class a showed higher expression levels at different stages of fiber and ovule development. The class b genes GhKNOX1-D, GhKNOX3b-A, GhKNOX6-A, GhKNOX10-A, and GhKNOX11-D showed higher expression in floral organs, such as the torus, sepal, and bract. Most class c genes showed higher expression levels in tissues except fibers. Among these genes, GhKNOX4b-D, GhKNOX12-D, and GhSTM2-D were more highly expressed in the root, whereas GhKNL1-A, GhKNOX4b-A, GhKNOX8-D, GhKNOX9-D, GhKNOX10-D, GhKNOX14-A/D, and GhSTM2-A were predominantly expressed in ovules. The class d genes showed diverse expression patterns, which were focused on the root, sepal, anther, and filament. These results indicated that GhKNOX genes may have diverse biological functions in different tissues.

The expression patterns of GhKNOX genes in the shoot meristem of the early maturing cultivar ‘Zao1’ and the late-maturing cultivar ‘CCRI50’ were analyzed from the fourth leaf expanded to the seventh leaf expanded stages (Supplementary Figure 3). Eight genes in class A showed decreased expression levels at the four shoot apical development stages of ‘Zao1’ compared with those of ‘CCRI50.’ The transcript level of class B genes was highest at the fourth leaf expanded stage of ‘Zao1’ and at the seventh leaf expanded stage of ‘CCRI50.’ Other GhKNOX genes in class C exhibited higher expression levels in ‘Zao1’ with a lower expression level detected at the fourth leaf expanded stage. In ‘CCRI50,’ the majority of GhKNOX genes showed the highest transcript level at the sixth and seventh leaf expanded stages except GhKNOX11-D. Six STM homolog genes showed higher expression levels in ‘Zao1’ than in ‘CCRI50.’ Thus, the functions of GhSTM genes in cotton growth and development require further verification.

The expression pattern of the 44 GhKNOX genes in response to exposure to cold, heat, salt, and drought stress was analyzed at different time points. The expression of some KNOX genes was affected significantly, such as GhKNL1-D, GhKNOX2-D, GhKNOX3b-A, GhKNOX4b-A, GhKNOX6-A, GhKNOX10-A, and GhKNOX14-A. The expression level of GhKNOX2-D, GhKNOX4b-A, GhKNOX6-A, and GhKNOX14-A was increased in response to the four stresses. The genes GhKNOX2-A, GhKNOX14-D, and GhKNL1-A showed decreased expression under the four stress treatments. Expression of GhKNOX10-A was not influenced by heat, drought, and salt stress. GhKNOX5a-D and GhKNOX7b-D showed higher expression levels in response to cold stress only, whereas expression of GhKNOX5a-D, GhKNOX9-A, and GhSTM3-D increased at 1 h and thereafter decreased slightly. The present results indicated that GhKNOX genes from the A subgenome displayed superior adaptability to environmental stresses (Figure 6).

To explore the expression of GhKNOX genes in response to abiotic stresses in greater detail, we selected eight GhKNOX genes for which expression was induced by drought and salt stress, and examined their expression following treatment with 20% PEG or 200 mM NaCl. The qRT-PCR results showed that GhKNOX4b-A/D, GhKNOX7b-A, GhKNOX10-A, and GhKNOX14-A were upregulated by PEG or NaCl treatment. Transcription of GhKNOX2-A and GhKNOX3b-A/D was upregulated by PEG treatment and downregulated by NaCl treatment. GhKNOX6-A, and GhKNOX9-A were downregulated by PEG or NaCl treatment (Figure 7). These results implied that GhKNOX family genes may show differential expression levels under different abiotic stresses.

GhKNOX2-A, GhKNOX10-A, GhKNOX14-A, GhSTM2-A, and GhSTM3-A/D belonged to the clade KNOXI, which includes the Arabidopsis homolog genes KNAT1, KNAT2, KNAT6, and STM. The expression patterns implied that GhKNOX2-A, GhKNOX10-A, and GhKNOX14-A are induced by salt stress. We used VIGS assays to investigate the functions of these five G. hirsutum genes. The appearance of white leaves indicated that VIGS was successful and qRT-PCR analysis confirmed that the expression levels of the five KNOX genes decreased significantly in the VIGS plants (Supplementary Figure 4). Silencing of GhKNOX2-A increased salt tolerance, therefore the silenced cotton seedlings grew better than the control seedlings in response to salt treatment (Figure 8A). The POD activity of the silenced plants was significantly higher than that for control seedlings (Figures 8G,H). Silencing of GhKNOX10-A and GhKNOX14-A decreased the salt tolerance (Figures 8B,C), therefore the silenced cotton seedlings showed inferior growth compared with the control seedlings in response to salt treatment. The MDA content in GhKNOX14-A VIGS plants was significantly higher than that in control seedlings, whereas the POD activity of silenced GhKNOX10-A plants was lower than that of control seedlings (Figures 8G,H). Compared with control plants, VIGS of GhSTM2-A and GhSTM3-A/D did not result in significant changes in MDA content after salt treatment, whereas the POD activity decreased compared with that of the control (Figures 8D–H). The flowering time was promoted in GhSTM3-A/D VIGS plants, and expression of GhFT and GhAP1 was upregulated with silencing of GhSTM3-A/D (Figures 8F,I). These results indicated that the five KNOX genes play an important role in salt stress tolerance and GhSTM3 might affect the floral transition of cotton.

The increase in availability of resequencing data for cultivated cotton species enabled assessment of genetic differences in KNOX genes over several decades of breeding. In this study, we estimated the genetic variation of 82 early and 67 modern cultivars that were sequenced and the data released from a core collection of upland cotton (Ma et al., 2018). The early cultivars included introductions and cultivars bred before 1976, and the modern cultivars comprised those bred during the period 1996–2008. To compare genetic variation among different KNOX family genes in cotton cultivars, we counted the number of SNPs per gene. A total of 64 SNPs were detected in 19 GhKNOX genes and the number SNPs per gene ranged from 1 to 11. The early cultivars contained 54 SNPs in 16 GhKNOX genes, whereas modern cultivars contained 57 SNPs in 18 GhKNOX genes. The SNPs density of modern cultivars was higher than that of early cultivars for the genes GhSTM1-A, GhSTM2-D, GhKNOX4b-D, GhKNOX8-A, and GhKNOX12-A, whereas the reverse result was observed for GhKNOX7b-D and GhSTM1-D. These results showed that the GhKNOX genes exhibited rich genetic variation among both early and modern cultivars (Figure 9A).

To clarify the selective pressure exerted during breeding, we estimated the genetic difference among the two groups of cultivars (Figure 9B). There were distinct selective signals for GhKNOX11-D (0.10), GhKNOX8-A (0.089), GhSTM1-D (0.65), GhKNOX2-A (0.51), and GhKNOX4a-D (0.48) during cotton improvement, showing that these genes were subjected to intensive artificial selection. In contrast, 11 genes (GhKNOX4b-A, GhKNL1-A, GhKNOX7a-A, GhKNOX2-D, GhKNOX7b-D, GhKNOX8-D, GhKNOX5a-D, GhSTM2-D, GhKNOX4b-D, GhKNOX7a-D, and GhKNOX1-D) showed few genetic differences and an average F_st of 0.006. These results indicated that the latter genes have not been subjected to breeding selection and are potential improvement targets for breeders in the future.