In order to identify members of the SNAT gene family, protein sequences and genomes of upland cotton were downloaded from the cotton database Gossypium Resource And Network database (http://grand.cricaas.com.cn) and CottonFGD(https://cottonfgd.org) (Zhu et al., 2017). The conserved domain of rice SNAT(AK059369) protein sequence was found to be Acetyltransf_1 by Pfam database (https://pfam.xfam.org/) analysis (Kang et al., 2013). Using the hidden Markov model of Acetyltransf_1 as a query file, genes with the conserved domain of Acetyltransf_1 in upland cotton genome were screened by the local software HMMER as candidate genes of SNAT gene family. Based on the conserved domain of SNAT protein, Acetyltransf_1(PF00583) was further analyzed by using the Pfam database, and genes with incomplete domain were manually deleted. Based on the location of the gene on the chromosome, we renamed the gene as GhSNAT1A-GhSNAT24A, GhSNAT1D-GhSNAT28D. The physicochemical properties of GhSNAT genes were analyzed using the online tool Expasy-Protparam (https://web.expasy.org/protparam/) (Gasteiger et al., 2005). In order to understand the subcellular localization of SNAT protein, we used online websites WOLF-PSORT (https://wolfpsort.hgc.jp/) for prediction.

In order to study the evolutionary relationship between SNAT genes, the protein sequence of upland cotton SNATs was input into software MEGA7, and ClustalW was used for multiple sequence alignment. The intraspecific evolutionary tree was constructed by neighborhood method, and the specific parameters were as follows: Bootstrap replication: 1,000, model/method: P-disrance, and all/Missing Data Treatment: Partial deletion.

In order to explore the evolutionary relationship of SNAT among different species, the online Phytozome V12.1 database was used to download the protein sequences of Theobroma cacao, Arabidopsis thaliana, Oryza sativa, Populus trichocarpa, Vitis vinifera, Glycine max, Zea mays. The homologous sequences of SNAT in these species were obtained by using the above methods. Protein sequences from four cotton species (Gossypium arboretum, Gossypium raimondii, Gossypium hirsutum, Gossypium barbadense) and members of the SNAT family of Theobroma cacao, Oryza sativa, Arabidopsis thaliana were input into MEGA7. ClustalW was used for multiple sequence alignment and Maximum Likelihood method was used to construct interspecific evolutionary tree. Parameter Settings were as follows: bootstrap replication 500, model/method: JTT + G, and all/Missing Data Treatment: Partial deletion (Kumar et al., 2016).

The gene structure of GhSNAT members was obtained by using the annotation file of whole genome of upland cotton downloaded from online cotton database CottonFGD. Online software MEME was predicted the motif of genes (http://meme-suite.org/tools/meme), and the parameters were as follows: the maximum number of motif of each gene is 20, and other parameters are set by default. The evolutionary relationship, gene structure and motif composition of SNAT were analyzed using TBtools (Chen et al., 2020).

The whole genome annotation files of four cotton species were downloaded from cotton database CottonFGD (https://cottonfgd.org/about/download/annotation), and gene replication events occurred in SNATs were analyzed by MCScanX software. The chromosome positions and gene replication of SNATs members of four cotton species were visualized by TBtools software (Chen et al., 2020).

The 2000 bp DNA sequence of the upstream region of GhSNAT was obtained from CottonFGD database (https://cottonfgd.org) as the promoter (Martinez and Chrispeels, 2003). We used PlantCARE (http://:/bioinformatics.psb.ugent.be/webtools/PlantCARE/html/) to predict cis-regulatory elements in the GhSNATs gene promoter region. Cis-acting elements related to plant hormones and abiotic stress were selected for further analysis and plotted by software TBtools (Chen et al., 2020).

In order to study the expression patterns of GhSNAT gene family, we obtained the expression levels of these genes under cold, heat, salt, PEG stress and the expression levels of these genes in roots, stems and leaves by Gossypium Resource And Network database online (http://grand.cricaas.com.cn/page/tools/expressionVisualization). Heat maps were generated using TBtool software using fragment number per kilobile exon (FPKM) (Chen et al., 2020).

In order to study the evolutionary relationship of SNATs in four cotton species, MCScanX software was used to analyze the collinearity of repeated gene pairs in Gossypium arboretum, Gossypium raimondii, Gossypium hirsutum, Gossypium barbadense (Wang et al., 2012), the software TBtools was used to visualize the results (Chen et al., 2020).

To determine selection pressure, Ka (non-synonymous substitution) and Ks(synonymous substitution) rates of duplicate genes were calculated using TBtools software (Chen et al., 2020).

Using upland cotton cultivar Zhong 9807 as experimental material, the seeds were sown on the medium of sand and vermiculite 1:1.5, and grew in an indoor incubator at 25°C for 16 h in the day and 8 h at night. In order to explore the effect of increasing melatonin level on GhSNATs expression, cotton seedlings at three-leaf stage were treated with 20 µM melatonin and sprayed on leaves once a day for three consecutive days. The seedlings were grown in an indoor incubator at 25°C for 16 h during the day and 8 h at night (Zhang et al., 2021a). Melatonin treatment with 0 µM concentration was used as control. The cotton seedlings treated with melatonin and the control group were treated with NaCl solution at 100 mM concentration for 12 h and sampled respectively. The three biological replicates were performed.

Total RNA was extracted by EASYspin Plus Plant RNA rapid isolation kit (Aidlab Co., LTD, Beijing, China). The pure RNA was reverse-transcribed using HiScript Ill RT SuperMix for qPCR (+gDNA wiper) (Vazyme Biotech Co., LTD, Nanjing, China) according to the manufacturer’s instructions. Gene-specific primers were designed using Primer Premier 5 software (Supplementary Table S1). qRT-PCR assays were performed on the Bio-Rad 7,500 fast fluorescence quantitative PCR platform with TransStart^® top green qPCR supermix (TransGene Biotech Co., LTD, Beijing, China) in accordance with the manufacturer’s protocol, three biological replicates, 2^−ΔΔCt methods were used to measure relative gene expression levels (Livak and Schmittgen, 2002). The internal control was GhUBQ7 (Tan et al., 2013), which was expressed stably in cotton plants and was not affected by treatment or genotype.

In order to investigate whether the content of endogenous melatonin in cotton affected the salt tolerance of cotton, the SNAT gene was silenced by virus-induced gene silencing. We used the published protein sequence of rice SNAT(AK059369) as the query sequence (Kang et al., 2013), the putative GhSNAT(GH_D02G1113.1) in cotton was obtained by comparing the cotton genome with local software BLAST, which was named GhSNAT3D in this study. VIGS vector was stored in our laboratory as pYL156 vector, BamHI and SacI restriction endonuclease were selected for double digestion, and online tool SGN-VIGS (https://vigs.solgenomics.net/) was used to design silent fragments with a length of about 300bp. The Infusion primer for GhSNAT was manually designed, and the primer for GhSNAT silent fragment was as follows: forward primer, 5′-TAG​AAG​GCC​TCC​ATG​GGG​ATC​CGA​ATT​TGT​GCT​TGT​TGA​AAA​GTC​TCA-3′ reverse primer, 5′-TGC​CCG​GGC​CTC​GAG​ACG​CGT​GAG​CTC​AGG​ATC​AAC​AAG​AAC​ATC​CCA​GA-3′. Using cDNA from cotton leaf tissue as template, the silent fragment was amplified and the VIGS expression vector pYL156:GhSNAT was constructed by in-fusion technique. The constructed expression vector was transformed into EScherichia coli. After correct sequencing, the expression vector was transformed into Agrobacterium tumefaciens by freeze-thaw method. The virus-mediated gene silencing (VIGS) system consisted of recombinant vector, negative control pYL156, positive control pYL156:PDS and helper vector pYL192. Therefore, the correctness of the whole system was judged by observing whether the true leaves of pYL156:PDS VIGS turned white. Prepare cotton seedlings for infection after cotyledons flatten, and water them thoroughly the night before infection. Make an incision in the lower epidermis of the cotyledon with a needle, and inject the bacteria into the lower epidermis of the cotyledon with a syringe until the bacteria fill the whole cotyledon. After the injection, the cotton seedlings were sheltered from light for 24 h and then cultured normally. Cotton was treated with 100 mmol L⁻¹ NaCl when it reached the three-leaf stage.

We took samples for the determination of endogenous Melatonin content, and measured the endogenous Melatonin content by Plant Melatonin (MT) ELISA Kit (Ziker, ZK-P7490, Shenzhen, China). This kit is intended for use only with a one-step sandwich ELISA kit. The 0.1 g sample was mashed with proper amount of normal saline, centrifuged at 3,000 rpm for 10 min, and the supernatant was taken. Melatonin was detected according to the Plant MT ELISA Kit, and each sample had three biological replicates.

Peroxidase (POD) activity and superoxide dismutase (SOD) activity were determined by using POD activity detection kit (Sinobestbio, YX-W-A502, Shanghai, China) and superoxide dismutase (SOD) activity assay kit (Sinobestbio, YX-W-A500-WST-8, Shanghai, China). About 0.1 g tissue was weighed and 1 ml extract was added for ice bath homogenization. Centrifuged 8,000 g at 4°C for 10min, supernatant was taken and placed on ice for determination according to the instructions, with three biological replicates for each sample.

The determination of Ca²⁺ content refers to EDTA titration method in GB 5009.922016 “Determination of Calcium in Food of National Standard for Food Safety”. Samples that needed to be assayed for Ca²⁺ content were cleaned by ddH₂O, placed in an oven, oven dried at 110°C for 10 min, 80°C until constant weight, accurately weigh 0.2 g of sample into a graduated digestion tube, add 10 ml 10% nitric acid, and digest in an adjustable electric furnace (reference conditions: 120°C/0.5 h to 120°C/1 h, increase to 180°C/2 h to 180°C/4 h, and increase to 200–220°C). The digestion solution appeared colorless and transparent or slightly yellow. Constant volume to 25 ml with water after cooling, then dilute as needed for the actual assay, and add a volume of lanthanum solution (20 µg/L) to the dilution to a final concentration of 1 μg/L and mix for further use, this being the sample to be tested. Pipette 1 ml of the sample to be tested and a blank into a test tube and add 1 drop of sodium sulfide solution (10 g/L), 0.1 ml of sodium citrate solution (0.05 mol/L), 1.5 ml of potassium hydroxide solution (1.25 mol/L), and 3 drops of calcium red indicator. Titrate immediately in a 10 fold dilution of EDTA solution until the indicator changes blue from purple red and record the volume of EDTA solution consumed with a 10 fold dilution.

In order to explore the interaction network of GhSNAT protein, protein sequences of GhSNAT3D were compared to Arabidopsis thaliana to obtain the homologous genes of GhSNAT3D in Arabidopsis. STRING database (https://string-db.org/) was used to construct the interaction network and analyze the functions of GhSNAT protein.

The GraphPad Prism 8.0 software was employed to analysis (ANOVA) the results. Duncan’s Multiple Range Test was used to compare the least signifcant difference of means (p < 0.05).

The latent Markov model of Acetylation transf_1 was used as a query file, and genes with conserved domain of Acetylation transf_1 were selected as candidate genes of SNAT gene family in upland cotton by HMMER software. Pfam database was used for further analysis, and genes with incomplete domain were removed manually. We renamed the genes GhSNAT1A-GhSNAT24A, GhSNAT1D-GhSNAT28D based on their location on the chromosome. We then analyzed and predicted the physical properties of these genes, including ID, isoelectric point, protein molecular weight, protein length and subcellular localization (Supplementary Table S2). A total of 52 SNAT genes were identified in the whole upland cotton genome, which encoded protein sequences ranging from 157 (GhSNAT1A, GhSNAT1D) to 610 (GhSNAT4A) amino acids. The isoelectric point ranged from 5.083 (GhSNAT9A, GhSNAT10D) to 10.08 (GhSNAT25D), and the molecular weight ranged from 17.718 (GhSNAT1D) kDa to 67.138 (GhSNAT16A) kDa. Subcellular localization predicted 21 genes in cytoplasm, 3 genes in cytoskeleton, 2 genes in mitochondria, 20 genes in chloroplast, 6 genes in nucleus.

In order to further study the evolutionary relationship of SNAT in upland cotton, we constructed an evolutionary tree of intraspects-specific SNAT genes in upland cotton (Figure 1A). According to the evolutionary relationship of genes and protein structure characteristics, the SNAT genes in upland cotton can be divided into five branches. Branch Ⅰ contains 24 genes and branch Ⅱ contains 13 GhSNATs, branch Ⅲ contains 4 GhSNATs, branch Ⅳ has 6 GhSNATs, and branch Ⅴ has 5 GhSNATs. In the phylogenetic tree, two genes gather together to form gene pairs, which form 22 gene pairs in total. In each pair, one gene comes from the A subgenome and one from the D subgenome.

To study the evolutionary relationships of SNAT in plants, 468 protein sequences from upland cotton, island cotton, Asiatic cotton, Raymond’s cotton, cocoa, Arabidopsis, rice, soybean, grape, poplar and maize (27 from Gossypium arboretum, 59 from Gossypium barbadense, 52 from Gossypium hirsutum, 26 Gossypium raimondii, 33 from Arabidopsis thaliana, 31 from Oryza sativa, 48 from Theobroma cacao, 49 from Zea mays, 102 from Populus trichocarpa, 114 from Glycine max and 36 from Vitis vinifera) were constructed to study the relationship between SNAT genes (Figure 1B). The SNAT proteins of these species were distributed in almost every clade. The phylogenetic tree was divided into five branches randomly distributed. Among these clades, clade Ⅲ had the least number of members (26), clade Ⅰ had the most members (141), clade Ⅱ, Ⅳ and Ⅴ contained 111, 85 and 70 genes, respectively. It should be noted that in Arabidopsis, rice and cacao, SNAT proteins of four cotton species all have corresponding homologous genes in clade Ⅰ - Ⅴ, indicating that the SNAT proteins of these plants are closely related in evolution. In phylogenetic trees, we found that the GhSNAT gene pairs and GbSNAT gene pairs were always clustered together, which could serve as evidence of gene duplication. At the same time, SNAT proteins of tetraploid cotton (Gossypium hirsutum, Gossypium barbadense) and diploid cotton (Gossypium arboretum, Gossypium raimondii) were clustered together, confirming that upland cotton and island cotton were the result of hybridization between Gossypium arboretum, and Gossypium raimondii.

To further understand the possible structural evolution of GhSNATs, we constructed a phylogenetic tree, gene structure, and motif association analysis of GhSNAT members (Figure 2). We constructed a phylogenetic tree using the protein sequences of GhSNAT members, obtained the gene structure of GhSNAT members from the whole genome annotation file of upland cotton, submitted the protein sequences of GhSNAT members to the online tool MEME for conservative motif prediction, and used TBtools for association analysis.

From gene structure and phylogenetic tree results, we can see that similar genes clustered together in the same set of phylogenetic trees. The number of exons in each gene ranged from 1 to 13. In most cases, two genes in the gene pairs had similar exon-intron structure and length, but some gene pairs were different in structure and length, such as GhSNAT4D/GhSNAT13D and GhSNAT4A/GhSNAT5D. There were 8 genes with only one exon and 2 genes with 13 exons. Different clades have different numbers of exons, but most GhSNAT members of the same clades share the same exon-intron structure. We found that GhSNAT gene structure is strongly related to phylogeny on an evolutionary basis. The motifs of each protein vary from 1 to 10. In addition, two genes in most gene pairs have the same motif composition, meaning that they are functionally similar at the protein level. And all GhSNAT proteins contained Motif1, a conserved domain shared by the GhSNAT family. Clade I mainly consists of motifs 1 and 2, and most members in clade I also contained motif5 and 11. Most of the GhSNAT motif compositions in clade Ⅱ were irregular, but all have motif1. The motif composition of clade III was more conserved and consists of motif5, 1, 15, 20, 4, 17. Clade Ⅳ mainly consists of motif12, 19, 3, 16, 8, 5, 6, 9, and 1, but GhSNAT22A and GhSNAT24D consist of motif4, 1, and 2; most members of clade V consist of motif10, 13, 4, 15, 7, 14, 1, 16. Interestingly, the GhSNAT13D/GhSNAT4D genes were not only different in gene structure, but also in motif composition, one of which was composed of Motif6, 1, 2, and one of which was composed of motif7, 1, 2. We speculated that it may have been mutated during evolution leading to changes in gene structure and function.

In order to study the chromosome distribution and gene replication of members of the SNAT gene family, we mapped the physical locations of these genes on cotton chromosomes (Figure 3). 164 SNATs were located on specific chromosomes of four cotton species. In Gossypium hirsutum, 52 genes were mapped onto 21 chromosomes, with uneven distribution. The number of SNAT genes on each chromosome was between 1-5. There were 24 genes in A subgenome and 28 genes in D subgenome. Chromosomes A03, A09, A12, D04 and D12 did not have SNAT. The number of genes on chromosomes A01/D01, A02/D02, A05/D05, A06/D06, A07/D07, A08/D08, A10/D10 was consistent, while the number of genes on other A/D chromosomes was different, which might be caused by the duplication or loss of SNAT members in the process of evolution. Interestingly, GhSNATs did not undergo tandem repetition during evolution, and fragment replication was the main mode of gene amplification.

In Gossypium barbadense, 59 genes were distributed on 22 chromosomes, with 1-6 SNATs on each chromosome, 29 SNATs on the A subgenome and 30 SNATs on the D subgenome. There were 6 SNATs on chromosomes A04 and D08, but no SNATs on chromosomes A09, A12, D04 and D12, suggesting that gene loss may be one of the reasons for species evolution. Although the number of SNAT in A/D subgenome was similar, the distribution of these genes in chromosomes was not corresponding. The number of genes on chromosome A01/D01, A05/D05, A06/D06, A07/D07, A10/D10 was the same, while the number of genes on other A/D chromosomes was different. In Gossypium barbadense, there were five tandem replication on chromosome A04 and one tandem replication on chromosome D02. Tandem replication was also one of the driving forces for the amplification of GbSNATs in evolution.

In Gossypium arboretum, 27 GaSNATs were unevenly distributed on 9 chromosomes, wherein chromosome Chr06, Chr09 and Chr10 have 1 GaSNAT respectively, and chromosome Chr08 and Chr11 have at most 5 GaSNATs. No genes were distributed on chromosome Chr01, Chr02, Chr04, and Chr12, and a tandem replication occurred on chromosome Chr03. In Gossypium raimondii, 26 GrSNATs were distributed unevenly on 10 chromosomes, and no GrSNATs were distributed on chromosomes Chr06, Chr08, and Chr12. In Gossypium raimondii, GrSNATs did not undergo tandem replication during evolution.

Promoters can interact with transcription factors to control the onset and degree of gene expression. Cis-acting elements are located in the promoter region of genes and can be used as a reference for stress response and tissue specificity in different environments. To understand the response mechanism of GhSNAT to abiotic stresses, we obtained FPKM of GhSNATs from Gossypium Resource And Network database online to analyze the expression patterns of these genes under various stresses (cold, heat, salt, and PEG). Phylogenetic tree analysis, promoter analysis and expression heat map association analysis were performed (Figure 4). Therefore, the analysis of GhSNAT promoter region was helpful to explore the potential function of genes. Cis-acting elements played an important role in stress response. We used 2000 bp DNA sequence of the upstream region of GhSNAT as the promoter, and predicted cis-acting elements of the promoter region of GhSNAT gene by using the online tool PlantCARE. A large number of cis-acting elements were detected in the promoter region, and cis-acting elements related to plant hormones and abiotic stresses were selected for further analysis.

In the plant hormone response category, the selected elements are abscisic acid response element, salicylic acid response element, gibberellin response element, MeJA response element and auxin response element; the elements responding to abiotic stress include defense and stress response element, wound response element, drought-inducibility element and low-temperature response element. More than half of GhSNAT members had MeJA-responsive, abscisic acid responsiveness, salicylic acid responsiveness, 41 genes contained abscisic acid responsiveness elements, 8 genes contained wound responsiveness, 5 genes contained gibberellin-responsive, and 13 genes contained auxin-responsive. There were 20 GhSNATs with defense and stress responsiveness elements, and there were 19 GhSNATs with drought-inducibility elements. 23 GhSNATs contained low-temperature responsive elements. We speculated that these genes are involved in abiotic stress responses along with hormo responses.

A total of 52 FPKM of GhSNATs were obtained from Gossypium Resource And Network database, and heat maps were made based on the expression levels of these genes under cold, heat, salt and PEG stresses. The results showed that only a few GhSNATs, such as GhSNAT14D, GhSNAT28D, GhSNAT23A, GhSNAT27D, GhSNAT20A were not differently expressed under various abiotic stresses. Most of the genes were strongly induced and differentially expressed under various stresses. We found that most of the genes from the same clade had the same expression pattern. Interestingly, some genes were induced by specific stresses, such as GhSNAT1D, which was strongly induced under cold stress but not under other stresses; GhSNAT9D was only induced by heat stress. We counted the number of differentially expressed GhSNAT genes under different stresses, including 29 differentially expressed GhSNATs under cold treatment, 35 differentially expressed GhSNAT genes under high temperature treatment, 40 differentially expressed GhSNAT genes under PEG treatment and 37 differentially expressed GhSNAT genes under salt stress. The expression levels of GhSNAT genes changed under different stress conditions, suggesting that GhSNAT members were involved in the regulation of abiotic stress.

Meanwhile, to explore the tissue specificity of GhSNATs, we obtained FPKM of GhSNATs expression in different tissues (root, stem, leaf, petal, torus) from Gossypium Resource And Network database And visualized the data by using TBtools. In the GhSNAT family, 7 genes (GhSNAT8D, GhSNAT21A, GhSNAT23D, GhSNAT28D, GhSNAT23A, GhSNAT18D, GhSNAT20A) showed low expression in different tissues. GhSNAT19A and GhSNAT21D were highly expressed in each tissue, and most of the genes were tissue specific, such as GhSNAT5D, GhSNAT11D, GhSNAT13A and GhSNAT2D were highly expressed in leaf. GhSNAT1D and GhSNAT2A were highly expressed in stem and petal. GhSNAT9D, GhSNAT13D and GhSNAT12A were highly expressed in torus.

Replication events, including whole genome replication, fragment replication and tandem replication, play an important role in gene amplification. Most plants have experienced ancient genome-wide replication events or polyploidy, and the duplication regions caused by genome-wide replication are usually the duplication of all genes on a large scale, rather than the duplication of a single gene or multiple genes (Malik et al., 2020). In biological evolution, the genome will undergo structural and quantitative changes during the processes of whole genome replication, chromosome recombination, chromosome inversion and translocation (Dujon et al., 2004; Nakatani et al., 2007). Large-scale genome-wide replication (WGD) and small-scale tandem replication and fragment replication between species can be identified from colinear fragments, which can be used as data for species tree inference. Through homology analysis of SNAT genes of four cotton species (Gossypium arboretum, Gossypium raimondii, Gossypium hirsutum, Gossypium barbadense), the loci relationships of SNAT genes of four cotton species were visualized (Figure 5). Genes on the same chromosome and adjacent to each other belong to tandem duplications, and the remaining genes from the same genome belong to segmental duplications (Liu et al., 2010), whereas the remaining genes from different genomes and subgenomes belong to the whole genome duplication. Gene duplication events are one of the major contributors to evolutionary dynamics and they have a major impact on genomic rearrangements and expansions (Chothia et al., 2003). We identified homologous gene pairs in GhA, GhD, GbA, GbD subgenomes of two tetraploid cotton and GaA, GrD genomes of two diploid cotton. Collinearity showed that there were several highly conserved loci in the A and D subgenomes of the two tetraploid cotton species, which was also the result of the cross between two tetraploids (Gossypium hirsutum, Gossypium barbadense) and two diploids (Gossypium arboretum, Gossypium raimondii).

Genes connected by lines of the same color represent the same gene. In Figure 5, we can see that the GhA/GhD, GbA/GbD subgenomes and many chromosomes in A and D genomes are connected by lines of the same color, that is, the GhA/GhD and GbA/GbD subgenomes have homologous genes of SNAT in A and D genomes. These results indicated that these genomes/subgenomes are evolutionarily related, and most SNAT genes are preserved in polyploid evolution. By comparing the genomes and subgenomes of Ga-Ga, Ga-Gb, Ga-Gh, Gb-Gb, Gb-Gr, Gb-Gh, Gr-Gh, Gr-Gr and Gh-Gh, a total of 608 lineal/parallel homologous gene pairs were identified, and 7 pairs of duplicated genes showed tandem duplication. There were 123 pairs of duplicated genes with fragment replication and the remaining 478 pairs of duplicated genes with whole genome replication. The number of Lineal/Parallel homologous gene pairs in Gh-Ga, Gh-Gr, Gh-Gb, Gb-Gr, and Gb-Ga replicated in the whole genome of SNAT gene was 70, 81, 165, 83, and 79, respectively. The colinear gene pairs of Ga-Ga, Gb-Gb, Gh-Gh and Gr-Gr were 8, 57, 52 and 6, respectively. Therefore, we hypothesized that genome-wide replication events and fragment replication events mainly occurred during the evolution of SNAT gene, leading to gene amplification.

In order to explore the effect of Darwin’s positive selection on the differentiation of duplicate genes, the present study analyzed the ratio of non-synonymous substitution rate (Ka) to synonymous substitution rate (Ks) in the coding region of all duplicated genes. In genetics, Ka/Ks represents the ratio between the non-synonymous replacement rate (Ka) and the synonymous replacement rate (Ks) of two protein-coding genes. This ratio determines whether there is selective pressure on the protein-coding gene. In general, there are pros and cons (and often cons) of natural selection because non-synonymous substitutions can cause amino acid changes that may alter the conformation and function of proteins and thus cause adaptive changes. TBtools was used to calculate the Ka/Ks ratio of these gene pairs (Figure 6), and the selection pressure of duplicate gene pairs could be inferred from the Ka/Ks ratio. The Ka/Ks values of 414 gene pairs were calculated in four cotton varieties. Generally, Ka/Ks = 1 is considered neutral selection, Ka/Ks < 1 is considered negative selection, and Ka/Ks > 1 is considered positive selection. There were 13 pairs of genes whose Ka/Ks values were all greater than 1, indicating that these genes had undergone positive selection in the evolution process. There were 340 pairs of genes with Ka/Ks values between 0 and 0.5, 61 pairs of genes with Ka/Ks values between 0.5 and 0.99, 96.86% of them with Ka/Ks values less than 1. This suggests that SNATs underwent intense purification and selection during evolution. In Ga-Gb, Gh-Ga, Gb-Gb, Gb-Gh, Gb-Gr and Gh-Gr, the number of gene pairs with Ka/Ks value greater than 1 is 2, 2, 2, 2, 3 and 2, respectively, indicating that these genes have been actively selected in the evolutionary process and have recently undergone rapid evolution. But whether it brings about harmful or beneficial traits remains to be studied.

Exogenous melatonin application has been reported to improve salt tolerance of cotton (Zhang et al., 2021a). In order to explore the effect of exogenous melatonin on GhSNATs expression, we detected the relative expression levels of melatonin on members of cotton GhSNAT clade Ⅱ family under salt stress (Figure 7). Under normal growth, some GhSNATs were not affected by exogenous melatonin, such as GhSNAT1D, GhSNAT3A, GhSNAT8D, GhSNAT9A, GhSNAT25D. The other genes in clade Ⅱ were down-regulated by exogenous melatonin. Under salt stress, GhSNAT1D, GhSNAT3A, GhSNAT6A, GhSNAT8D, GhSNAT9A, GhSNAT25D were not affected by exogenous melatonin. GhSNAT1A and GhSNAT3D were down-regulated by melatonin. GhSNAT7A, GhSNAT7D, GhSNAT10D, GhSNAT23D, and GhSNAT25D were up-regulated.

In order to study the function of upland cotton Serotonin N-acetyltransferase (GhSNAT), melatonin synthesis gene GhSNAT3D was silenced by VIGS technique. PYL156:PDS can change the leaf from green to white. The success of the experiment can be judged by observing whether the leaves of pYL156:PDS transformed plants turn white or not, and the no-load plants of pYL156 were injected as blank control. At about 2 weeks after VIGS infection, the leaves of pYL156:PDS plants became albino, indicating that our silencing system was stable. The expression level of GhSNAT3D was detected at the three-leaf stage of cotton seedlings, and it was found that the expression level of GhSNAT3D was significantly decreased in pYL156 as a control (Figure 8A), indicating the successful gene silencing.

PYL156 and pYL156:GhSNAT3D were then subjected to salt stress using 100 mM/L NaCl solution, with pYL156 as control (CK) and pYL156:GhSNAT3D as V-S. After 3 days of salt treatment, the plants showed a distinct phenotype (Figure 8B). Compared with CK, V-S plants were more seriously stressed, and the true leaf wilting degree of V-S plants was significantly higher than that of control plants. Because GhSNAT3D is a putative melatonin synthesis gene, endogenous melatonin levels were significantly reduced in silent plants (Figure 8C). It has been reported that melatonin enhances salt tolerance of cotton mainly by affecting reactive oxygen species scavenging system and Ca²⁺ signal transduction pathway (Zhang et al., 2021a). POD activity, SOD activity and Ca²⁺ content of biochemical indexes were measured on plants CK and V-S respectively (Figure 8C). The POD activity, SOD activity and Ca²⁺ content of GhSNAT3D gene silenced plants were significantly down-regulated by gene silencing. We detected the relative expression of some Ca²⁺ signal transduction-related genes (Figure 8D). Compared with control plants, most Ca²⁺ signal transduction-related genes were up-regulated in GhSNAT3D gene silenced plants, as shown in: CAM, CAMBP25, CBL4, CBP60B, CBP60B1, CML8, CML10, CPK2, CPK5; GhSNAT3D gene silencing up-regulated the expression of CAL1 gene. CML7 was not induced to be differentially expressed. These results indicated that GhSNAT3D silenced plant activated oxygen scavenging system weakened, Ca²⁺ content decreased, Ca²⁺ signal transduction related genes were differentially expressed, and salt tolerance of cotton decreased.

GhSNAT silenced plants decreased endogenous melatonin content and increased salt sensitivity. In order to investigate whether exogenous melatonin supplementation can improve the effect of GhSNAT3D gene silencing on plants, we treated part of pYL156:GhSNAT3D plants with exogenous 20 μM/L melatonin and the same treatment. PYL156 plants were used as control (CK), pYL156:GhSNAT3D plants as V-S, pYL156:GhSNAT3D plants + melatonin as V-S'. After 3 days of salt treatment, V-S plants were more severely stressed than V-S ′plants, with cotyledons and true leaves wilting significantly in V-S plants, while cotyledons in V-S′ plants wilted slightly and true leaves lost luster (Figure 8). We measured POD activity, SOD activity and Ca²⁺ content of biochemical indexes, and found that exogenous melatonin increased POD activity, SOD activity and Ca²⁺ content of V-S plants (Figure 8). Exogenous melatonin supplementation alleviated the effect of GhSNAT3D gene silencing on salt tolerance of cotton, and alleviated the salt sensitivity of GhSNAT3D gene silencing plants.

To further understand the function of GhSNAT protein, we compared GhSNAT3D protein to Arabidopsis thaliana, and obtained the homologous gene SNAT/NSI of Arabidopsis thaliana (AT1G32070.1) (Lee and Back, 2018), used the online tool STRING to predict the interaction protein of SNAT protein (Figure 9). Among these interacting proteins, AT4G35160, AAS (AT2G20340) and AT1G26220 have been identified in curated Databases to interact with AtSNAT. 19 proteins were predicted to interact with AtSNAT by textmining. Seven proteins may interact with AtSNAT through co-expression analysis. There was a close relationship between AT1G26220 protein and AtSNAT in the following aspects: co-expression, curated databases, textmining and cooccurence. Interestingly, the protein AT1G26220 has been identified as AtSNAT2 in Arabidopsis thaliana (Lee et al., 2019), a member of the SANT family, and its homologous gene in land cotton is GhSNAT25D. Based on these results, we speculated that GhSNAT3D and GhSNAT25D regulated melatonin levels in cotton through their interactions. Among these interacting proteins, the proteins in the Melatonin BioSynthetic Process (GO:0030187) were AT1G26220 and AT4G35160. AT4G35160 has been identified as AtASMT, the last gene for melatonin synthesis in Arabidopsis thaliana (Byeon et al., 2016). At the same time, the interacting proteins AT1G26220, AT4G35160 and AAS were enriched in Tryptophan Metabolism (ATH00380), and the substrate for the synthesis of melatonin was Tryptophan. Therefore, GhSNAT3D regulates the level of cotton melatonin through a complex pathway, that is, GhSNAT3D interacts with GhSNAT25D, ASMT and AAS to participate in the synthesis of melatonin.