Transcriptomic profiling of cotton leaves in response to cotton aphid damage

Cotton aphid (Aphis gossypii Glover) is regarded as one of the most harmful insect pests for cotton production. They are usually capable of causing severe yield loss through sucking plant sap, secreting honeydew and transmitting plant viral diseases. However, the molecular mechanism of the interaction between cotton and cotton aphids remains unclear currently. Therefore, the RNA-Seq study of cotton leaves was performed in response to cotton aphid damage at different time points (0 h, 6 h, 12 h, 24 h, 48 h and 72 h). A total of 9, 103 new genes were identified, and 7, 510 of them were annotated functionally. Based on the comparison results, the gene expression was analyzed according to the expression amount of genes in different samples. 24,793 differentially expressed genes (DEGs) were authenticated in all and their functional annotation and enrichment analysis were conducted. Compared with 0 h (without aphid damage, CK), the amount of down-regulated DEGs was largely more than that of the up-regulated DEGs at different time points under cotton aphid attack except for 48 h. As revealed by the functional annotation of DEGs, these genes are involved in all kinds of plant biological processes, including various resistance to abiotic and biotic stress, hormone metabolism, signaling transduction and transcriptional regulation. The results revealed the molecular mechanism of the interaction between cotton and cotton aphids and would facilitate the development of plant aphid-resistant cultivars.


Introduction
Aphids are one of the common pests on almost all crops. They feed on plant phloem sap with their piercing-sucking mouthparts, which can have a direct impact on the normal growth and development of crops and eventually cause a severe loss of crop yield and quality. Compared with chewing insects, aphid feeding only causes relatively little mechanical injury to crops. However, the everlasting feeding of a large quantity of aphids usually consumes a large amount of plant resources. Secondly, a lot of aphid species have the potential to spread plant viruses. For them, the detrimental effects caused by virus transmission are likely to outweigh the direct effects of aphid feeding. In addition, sooty mould grows on honeydew secreted by aphids on the leaf surface to constrain leaf photosynthesis (Züst and Agrawal 2016;Ng and Perry 2004). To counter Communicated by Z.-L. Zhang.
Xiao Zhong and Yazhen Yang contributed equally to this work and should be considered as co-first authors. these damages, plants have developed a series of intricate mechanisms to combat against aphids. Plenty of studies have indicated that various changes occur in plants damaged by aphids, such as protein phosphorylation, calcium flux, reactive oxygen species (ROS) generation and phytohormone changes, which result in relevant transcriptional regulation in the early response to aphids. Finally, plants are capable of producing various defense compounds, including nutrient compounds, glutathione S transferases (GSTs), peroxidases, and secondary metabolites for the sake of self-defense (Douglas 2018). Overall, plant susceptibility or resistance to aphids is reliant on their ability to recognize aphid feeding and rapidly initiate defense response . Plant resistance is considered as one of the most effective ways of achieving plant health management at present and in the near future (Li et al. 2019). Therefore, constantly exploring new aphid-resistant genes has been a long-standing hot topic for scientists. Mi-1.2 from the tomato was discovered to contribute to resistance to Macrosiphum euphorbiae and Bemisia tabaci (Rossi et al. 1998;Nombela et al. 2003). The melon gene Vat was reported to safeguard against virus aphid transmission (Aphis gossypii Glover) (Villada et al. 2009). Both of them are the earlier known aphid resistance genes and fall into NBS-LRR family, indicating that NBS-LRR family members perform a significant part in plant aphid resistance. Lectins have also been proven to perform resistance to aphids. For instance, the Galanthus nivalis agglutinin (GNA) gene encoding a monocot mannose-binding was extensively and clearly documented to confer resistance to Myzus persicae (Hilder et al. 1995). Amaranthus caudatus agglutinin (ACA) gene conferred resistance to Aphis gossypii and Myzus persicae Guo et al. 2004). In addition, the other aphid-resistant genes, including rag, chitinase, ɑ-amylase and on the likes, are researched actively and widely (Yates and Michel 2018).
RNA-Seq is a robust technology for analyzing transcriptome and excavating new functional genes. Therefore, it is used in many fields including plant defense response to phytophagous insects (Rai et al. 2018). For instance, the changes in gene expression in the whole genome of the cucumber aphid-resistant cultivar 'EP6392' were monitored using an Illumina Genome Analyzer platform (Liang et al. 2015). Furthermore, a comprehensive insect resistance response mechanism in cotton infested by the phloem-feeding insect Bemisia tabaci was discovered through transcriptome analysis (Li et al. 2016). In addition, the transcriptome analyses of Gossypium hirsutum at 24 h and 48 h in response to aphid and whitefly were conducted comparatively (Dubey et al. 2013). However, compared to other research areas, RNA-Seq analyses of plant defense response to phytophagous insects have been relatively sporadic up to now.
As an essential economic crop, the recent research on cotton transcriptome characterization is focused on comparative transcriptomic analysis at different growth stages and under all sorts of stresses. For instance, many new genes of cotton fiber were extensively analyzed from different angles using RNA-Seq method (Yoo and Wendel 2014;Naoumkina et al. 2015;Yang et al. 2006;Liu et al. 2016;Ma et al. 2016). Furthermore, cotton RNA-Seq data were also analyzed under various abiotic stress, such as salt stress (Peng et al. 2014;Shu et al. 2017;Zhang et al. 2016;Zhu et al. 2018), dry stress (Bowman et al. 2013;Fracasso et al. 2016), and high-temperature stress (Min et al. 2014). Cotton is also attacked by all kinds of biotic factors, such as cotton Verticillium wilt. As one of the most destructive cotton diseases, the study on cotton Verticillium wilt using RNA-Seq is concerned widely. For example, RNA-Seq transcriptional analysis and histo-chemistry have illustrated that lignin metabolism plays an important role in the cotton resistance to Verticillium dahlia (Xu et al. 2011). The gene families of resistance gene analogues in cotton and their response to Verticillium wilt were also studied . The analysis of sea-island cotton and upland cotton in response to the damage of Verticillium dahliae was conducted using RNA-Seq (Sun et al. 2013). However, RNA-Seq has rarely been applied to research the interaction between cotton and herbivorous insects.
It is well known that using cotton aphid resistance to control aphids is more effective and eco-friendly. Therefore, the analysis of the transcriptome dynamics of cotton defense response to cotton aphid damage was conducted using the RNA-Seq method for a better understanding of cotton aphid resistance mechanism. Our research will lay a solid foundation for further research into the molecular mechanism of the interaction between cotton and aphids, searching for new cotton aphid-resistant genes and developing new cotton aphid-resistant varieties.

Plant materials and growth conditions
Red-leaf cotton (Gossypium hirsutum), an aphid-resistant cotton cultivar, was provided by the Institute of Cotton from the Chinese Academy of Agricultural Sciences. Cotton seeds were sown directly in 30-cm-diameter plastic pots filled with nutrient soil until the seedling grew to 4 leaves in a growth room (25 ℃, 16 h light/8 h dark). 30 adults of cotton aphids (Aphis gossypii) were transferred to the top leaves with a small soft brush, and the double-sided adhesive tape was used to wind the petioles of the top leaves to prevent aphids from escaping. The top leaves were collected, respectively, without cotton aphid damage (0 h, CK) and different plants, respectively, at 6, 12, 24, 48 and 72 h under the attack of cotton aphids, transferred quickly into liquid nitrogen for 1 min, and stored in the ultra-low temperature freezer.

Library construction and sequencing
Total RNA from each sample was isolated using the Trizol method. 1.5 μg RNA from each sample was used as input material for the removal of rRNA utilizing the Ribo-Zero rRNA Removal Kit (Epicentre, Madison, WI, USA). The libraries for sequencing were established with NEBNext R Ultra™ Directional RNA Library Prep Kit for Illumina R (NEB, Ipswich, USA). The index codes were appended to attribute sequences to each sample. The quality of the library was estimated using the Agilent Bioanalyzer 2100 and qPCR method. The clustering of each index-coded sample was executed on acBot Cluster Generation System using TruSeq PE Cluster Kitv3-cBot-HS (Illumia). After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform and 150-200 paired-end reads were derived.

Transcriptome assembly
Raw reads were pre-processed for quality filtering where reads including that (1) unknown (N) bases are more than 10%, (2) contain adaptor sequences, and (3) contain lowquality bases (Q ≤ 20) more than 50% were removed through in-house perl scripts. According to the reference genome originated from the CottonGen database (http:// www. cotto ngen. org), the mapped reads were assembled using Cufflinks software. The original unannotated transcription area was searched using BLAST software that was used to conduct sequence alignment with databases of NR (Leng et al. 2013), swiss-prot (Deng et al. 2006), GO (Apweiler et al. 2004), COG (Ashburner et al. 2000) and KEGG (Tatusov et al. 2000) and all parameters are default values, and new transcripts and new genes were obtained to supplement and perfect the original genome annotation information. The sequences that are less than 150 bp or that contain only a single exon were filtered out. The sequencing data were deposited in the NCBI Short Read Archive database with the accession number PRJNA576973.

Differentially expressed gene analysis
Gene FPKMs (fragments per kilo-base of exon per million fragments) were calculated by summing the FPKMs of transcripts in each gene group. FPKMs of coding genes in each sample were computed using StringTie (1.3.1). EBSeq was employed for different analysis (Kanehisa et al. 2004). During the course of differentially expressed gene detection, Fold change ≥ 2 and False Discovery Rate (FDR) < 0.05 served as the screening criteria. Fold change represents the ratio of expression quantity between two treats. FDR is determined through the p-value correction of different significance.

Functional annotations of differentially expressed genes
The functions of differentially expressed genes were annotated with NR, swiss-prot, GO, COG and KEGG. To obtain the annotation information about the new genes, the reference Gossypium hirsutum genome and the annotation files were downloaded from the CottonGen database (http:// www. cotto ngen. org). Gene Ontology (GO) enrichment analysis of DEGs was implemented by the topGO R packages. KOBAS software was used to test the statistical enrichment of DEGs in KEGG pathways.

Quantitative real-time PCR
To ascertain the expression levels of ten genes, quantitative real-time PCR (qPCR) was carried out based on the manufacturer's guides for the Bio-Rad CFX 96 Real-Time Detection System (Bio-Rad, Hercules, CA, USA) and the SYBR premix ex Taq II system (Takara perfect real time). Cotton leaves were firstly collected without cotton aphid damage (0 h, CK) and at 6, 12, 24, 48, and 72 h under cotton aphid attack at the four-leaf stage. Total RNA was isolated from each sample and reversely transcribed into cDNA as qPCR template. The expression profiles of the relative genes were examined through qPCR method. All reactions were conducted in triplicate, and controls were included. The 2 −△△Ct method was applied to the calculation of relative gene expression values (Pfaffl 2001). The prime sequence was presented in Table S1. The difference between each treatment (6, 12, 24, 48 or 72 h) and the control (0 h) was statistically analyzed with the t test for independent samples for qPCR data. Two significance levels were used (*, P < 0.05 and **, P < 0.01).

RNA-Seq outline of cotton leaves in response to cotton aphid attack
The damage caused by cotton aphids is rated as serious at the seedling stage during the growth and development of cotton. Therefore, the cotton leaves at the four-leaf stage were selected as the research materials. To excavate the crucial genes at the transcriptional level in cotton defense response to cotton aphid damage, the whole-transcriptome RNA sequencings of cotton leaves at 0 (without aphid damage, CK), 6, 12, 24, 48 and 72 h under sustained cotton aphid attack were performed, respectively. After filtering out reads containing adapter or ploy-N, and low-quality sequences, there were 60. 05, 58.70, 57.74, 65.73, 59.80, and 60.26 million clean reads with pair-end (Q 30 > 85.03%, GC% = 41.3-42.9%), respectively, produced by RNA-Seq in all the six libraries (Table S2). For the clean reads according to the single-end, 89.94-90.78% of them were mapped to the cotton reference genome sequences (Table S3). The mapped reads were assembled with Cufflinks software, and compared with the original genome annotation information to explore the original unannotated transcription area and excavate the new transcripts and new genes. A total of 79,581 unigenes were identified by filtering out sequences that were comprised of less than 150 bp or only a single exon. To obtain the annotation information of the new genes, the databases were analyzed using a combination of BLAST (Altschul et al. 2017), NR, GO, KEGG, COG and Swiss-Prot. A total of 9,103 new genes were discovered, of which 7,510 were functionally annotated (Table 1). The normalized FPKM was used to quantify the gene expression level (Trapnell et al. 2010). The transcript expression level (FPKM) in this study ranged from 10 -2 to 10 4 , which was found to be coherent with most of the RNA-Seq results (Fig. 1a). FPKM box plot analysis demonstrated that the gene transcript expression levels varied in the six RNA-Seq data to some extent (Fig. 1b).

Differentially expressed genes from cotton leaves under cotton aphid attack
As gene expression has time specificity, it is quite significant to study the cotton differentially expressed genes (DEGs) at different time points under cotton aphid attack. Each sample damaged by cotton aphids was compared with the control (0 h) to identify differentially expressed genes (6 h vs. control, 12 h vs. control, 24 h vs. control, 48 h vs. control and 72 h vs. control). Fold change ≥ 2 and FDR < 0.05 were taken as the screening criteria. From these pairwise comparisons drawn with the control sample, we identified 5,790 (2,580 up-and 3,210 down-regulated), 9,726 (3,452up-and 6,274 down-regulated), 5,279 (1,836up-and 3,443 down-regulated), 4,549 (2,477 up-and 2,072 downregulated) and 8,683(2,769 up-and 5,914 down-regulated)  (Fig. 2a, Table S4). The differentially expressed genes were analyzed using hierarchical clustering, and the genes with the same or similar expression patterns were clustered (Fig. 2b). In addition, it was discovered that the DEGs between cotton aphid attack (48 h) and CK (0 h) is distinct from the other treats. The amount of up-regulated genes was higher than that of down-regulated genes between cotton aphid attack (48 h) and CK (0 h) (Fig. 3a, b).
In conclusion, these results suggested that cotton's defense response to the damage of cotton aphids involved multiple genes and metabolic pathways, and cotton plants end up showing a comprehensive defense potential against the attack of cotton aphids.

Differentially expressed transcription factors involved in response to cotton aphid stress
Transcription factors (TFs) are specific DNA-binding protein that has regulatory functions at the transcription level and performs an essential role in plant growth and development. In this study, a total of 945 differentially expressed transcription factors were identified from cotton The parenthesis after each term label is the number of genes containing significant differences for that term. b The KEGG rich cluster class diagram of the total differentially expressed genes. Red indicates high expression functional classification, while blue indicates relatively low expression metabolic pathways. The color value is log2(FPKM + 1). Each metabolic pathway is labeled with the number of genes containing significant differences within the metabolic pathway in parentheses leaf response to cotton aphids. They were classed into 27 TF families, respectively (Table S6). A large majority of them encoded the members of the bHLH, MYB, WRKY, ERF, bZIP, TCP and NAC TF families (Fig. 5a). The bHLH family with 199 DEGs was rated as the largest TF family in cotton leaf response to cotton aphid stress. In addition, a total of 152 DEGs were identified to belong to the MYB family. The expression levels for many of the differentially expressed TFs were up-regulated under cotton aphid attack, especially at 48 and 72 h under cotton aphid attack. However, the transcripts of some differentially expressed TFs were down-regulated, indicating that the functions carried out by cotton TFs were complex in cotton response to cotton aphid attack (Fig. 5b).

Differentially expressed aphid-resistant genes involved in response to cotton aphid stress
Aphid damage often triggers the differentially expressed genes related to plant resistance, involving plant morphological resistance, phytoprotective enzyme and SA, JA and ethylene signal pathway. Plenty of cotton aphid resistance genes were identified at first in this study. For example, 39 differentially expressed genes of the NBS-LRR family were identified and 10 genes of them were reported at first. They were associated closely with the detection of bacteria, viruses, fungi, nematodes, insects and oomycetes (Table S7). 106 DEGs of leucine-rich repeat receptor-like protein kinase were also found and 6 genes of them were novel (Table S8). Besides, 10 genes of 140 DEGs related to lectins were regarded as cotton new genes. Compared to the expression levels without cotton aphid attack (0 h), the expression levels of most lectin genes were up-regulated at 6, 12, 24, 48 and 72 h, and only transcription of two lectin genes was down-regulated ( Fig. 6a, b, Table S9). Additionally, in the functional analyses of differentially expressed genes, we also obtained many new other genes related to cotton aphid resistance, including 3 genes related to callose (Table S10), 3 genes related to trichomes (Table S11), 3 genes related to waxes (Table S12), 7 genes related to ca 2+ (Table S13), 1 gene related to peroxidase (Table S14) and 1 related to Phenylalanine ammonialyase (Table S15). 1 NPR1 gene and 2 NPR3 genes related to SA signaling pathway, 13 lipoxygenase genes related to JA signaling pathway, and 101  (Table S16).

Verification of RNA-Seq data by qPCR
To validate the results of the RNA-Seq, ten random DEGs including five up-regulated and five down-regulated genes were selected to perform qPCR. The selected genes were successfully amplified and the products were of the expected size, indicating the reliability of the assembly work. The qPCR results demonstrated that the relative expression levels of ten selected genes were consistent with the results of RNA-Seq despite the differential expression folds, suggesting that the RNA-Seq data were highly reliable (Fig. 7).

Discussion
RNA-Seq presents an effective way for the identification of DEGs and their regulatory mechanisms at the transcriptome level (Wang et al. 2009). It could provide new insights into the molecular basis for plant response to all kinds of abiotic and biotic stress (Sjokvist et al. 2019;Ashoub et al. 2018;Janiak et al. 2018).

RNA-Seq in cotton response to cotton aphid damage
Aphids are the largest group of phloem-feeding insects and often cause severe economic loss in crop production. However, the mechanism for the interplay between plants and aphids remains unclear at the moment (Züst and Agrawal 2016). Therefore, there is great significance for understanding their interaction mechanism to study the cotton response to cotton aphid damage using RNA-Seq method. Dubey et al. (2013) carried out comparative transcriptome analysis of Gossypium hirsutum in response to aphids and whiteflies (Dubey et al. 2013). However, only the transcriptomes of cotton infested with aphids and whiteflies for 2 h and 24 h were sequenced and analyzed. In this study, the RNA-Seq data of cotton leaves at 0, 6, 12, 24, 48, and 72 h under the attack of cotton aphids were analyzed. The sequencing results suggested that a total of 9, 103 new genes were discovered and 7, 510 were annotated functionally at different time points. Based on the comparison results, the gene The number of downregulated DEGs was found out to be largely higher than that of the upregulated genes at a different time under cotton aphid attack, which was in line with the previous studies that infestations of whiteflies and aphids drove transcriptional suppression over induction (Dubey et al. 2013;Vos et al. 2005). However, the number of down-regulated DEGs was slightly smaller than that of the upregulated genes at 48 h under cotton aphid attack in this study. Functional annotation of DEGs led to a finding that these genes involved in a lot of plant biological processes, including all sorts of resistance to abiotic and biotic stress, hormone metabolism, signaling transduction and transcriptional regulation. However, only two genes in the "response to insect (GO: 0,009,625)" and the "defense response to insect (GO: 0,002,213)" were enriched in the biological process, which suggested that the study of biological information on plant response to insects should be enforced. In addition, the study on cotton transcriptome in response to the damage of aphids and whiteflies displayed significant enrichment of the amino acid biosynthesis pathway (Dubey et al. 2013). We got similar results in our study, suggesting the amino acid biosynthesis pathway was important very much.

Plant aphid-resistant genes
Plant defenses against aphids are initiated at various levels of their interplay with aphids (Douglas 1998). Firstly, the plant surface is the first barrier of aphid feeding and activities, such as thorn and glandular trichomes on the surface of the plant (Wójcicka 2015;Stoner 1990;Ellis et al. 1996). In our study, 3 new DEGs (Protein WAX2 or Protein WAX2 -like protein) were identified, and these DEG expressed levels varied at 6, 12, 24, 48 and 72 h after cotton aphid attack in comparison with the CK Control (0 h). Trichomes on the surface of plants could also exert impacts on aphid activities by hindering aphid movement. Furthermore, glandular trichomes are usually the origin of sugar esters and secondary metabolites that are harmful to insects (Glas et al. 2012;Goffreda et al. 1990). Some glandular trichomes even have the potential to release the aphid alarm pheromone, such as (E)-β-farnesene (Gibson and Pickett 1990). Our data analysis showed that a total of 3 differentially expressed genes on trichomes were identified at first, suggesting that trichomes were crucial to plant aphid resistance. In addition, plants inhibit the afflux of phloem inclusion into the aphid by facilitating sieving element occlusion (SEO). SEO includes two processes: a rapid formation of proteinaceous plugs that transiently seal sieve plates and a slower callose deposition that causes long-term occlusion of sieve tubes (Verma and Hong 2001). A total of 3 new DEGs on callose deposition were identified in this study. These results indicated that callose involved in cotton response to cotton aphid attack. Some loci with specific aphid resistance have been stated in many plants, but only the Mi-1.2 gene from tomato and the Vat gene from melon have been cloned. They contribute to resistance against Macrosiphum euphorbiae and Aphis gossypii, respectively (Rossi et al. 1998;Vos et al. 1998;Dogimont et al. 2014). They belong to coiled-coil nucleotide-binding site -leucine-rich repeat type proteins (CC-NBS-LRR). In our study, 10 new NBS-LRR DEGs were found and 6 Leucine-rich repeat receptor-like protein kinase were identified at first. The carbohydrate-binding protein lectins exist in many plants and perform important functions in guarding plants against insect pests. At the moment, many lectin genes from plants have been transferred to plants to improve resistance to many phloem insects (Duan et al. 2018). The agglutinin gene is the most widely studied for plant resistance to aphids as carbohydrate-binding proteins currently.10 genes were considered as cotton new lectin genes in our RNA-Seq data. In addition, Plenty of other aphid-resistant genes also involved in this study, including peroxidase, phenylalanine ammonialyase, polyphenol oxidase and the likes.

Transcription factors involved in response to cotton aphid stress
Some transcription factors in plants have been validated to act as a key role in the face of attack by herbivorous insects. For example, overexpressing CmMYB19 improved aphid tolerance in Chrysanthemum by facilitating lignin production (Wang et al. 2017). Similarly, three MYB genes were related to the wheat defense against English grain aphid (Zhai et al. 2017). In addition, overexpression of GsMYB15 from the wild Soybean R2R3 enhanced resistance to Helicoverpa Armigera in Arabidopsis (Shen et al. 2018). However, myb102 from Arabidopsis increased plant susceptibility to aphids by substantial activation of ethylene biosynthesis (Zhu et al. 2018). A total of 152 DEGs belonging to the MYB family were identified in this study and their functions need to further be verified through relative experiments. Some members of the WRKY family are reported for aphid resistance. For instance, the overexpression of a chrysanthemum WRKY transcription factor enhanced aphid resistance . However, AtWRKY22 from Arabidopsis increased Arabidopsis susceptibility to aphids and modulated SA and JA signaling (Kloth et al. 2016). In addition, the rice WRKY53 suppressed herbivore-induced defenses by acting as a negative feedback modulator of map kinase activity . In short, the functions of transcription factors in plant response to aphids are relatively complex. Studies of cotton transcription factors are primarily focused on fiber development and the response to all kinds of abiotic stress at present. Therefore, studies of their aphid resistance should be performed continuously in the future.

Plant defense signaling in cotton response to cotton aphids
The plant hormones (salicylic acid, jasmonic acid, ethylene and abscisic acid) lead to signaling related to plantaphid interaction. Despite the aphid effectors remain to be identified, important advances have been acquired in comprehending the signaling machinery related to Mi-1.2 in tomato. Some genes in ETI to microbes are also required for Mi-1.2-conferred resistance to Macrosiphum euphorbiae. Such as Heat shock protein 90 and Suppressor of G-two allele of Skp1 (Bhattarai et al. 2007). A receptor-like kinase encoded by the tomato SERK1, a mitogen-activated protein kinase cascade and the transcription factors WRKY70 and WRKY72 are also required for Mi-1.2 resistance to aphids (Atamian et al. 2012;Bhattarai et al. 2010;Mantelin et al. 2011). Ca 2+ functions as a secondary messenger in eukaryotes. It is universally thought that Ca 2+ could influence callose deposition, promote phloem protein aggregation and result in phloem occlusion (Vincent et al. 2017;Aidemark et al. 2009;Furch et al. 2009). 7 new DEGs on Ca 2+ were identified in this study, suggesting that ca 2+ was essential in cotton response to cotton aphids. However, aphids conquered this defense by secreting Ca 2+ -binding protein, which suggested that the interaction between plants and aphids was relatively complicated (Dubey et al. 2013).

Conclusions
In the present study, RNA-Seq was applied to the detection of the global transcriptional changes in seedling cotton leaves in response to cotton aphid attack. A total of 9, 103 new genes were identified and 7, 510 were annotated functionally. Crucially, 24,793 DEGs were commonly identified in response to cotton aphids, suggesting that there were plenty of common and unique molecular mechanisms in relation to the response to cotton aphid stress cotton. This study extends the understanding of the molecular mechanisms for cotton leaf resistance to cotton aphids in the seedling stage and facilitates the development of plant aphidresistant cultivars.
Author contribution statement ZJM and YYZ conceived and designed the research. ZX, FP, MQQ performed the experiments, ZJM analyzed the data and wrote the manuscript. SQ and WXP revised the manuscript. All the authors have read and approved this current manuscript.