Huan Si,Hongling Liu,Yiwen Sun,Zhongping Xu,Sijia Liang,Bo Li,Xiao Ding,Jianying Li,Qiongqiong Wang, Lin Sun, Xianlong Zhang, Shuangxia Jin

National Key Laboratory of Crop Genetic Improvement,Huazhong Agricultural University,Wuhan 430070,Hubei, China

Keywords:Oral secretions Helicoverpa armigera Spodoptera litura Cotton Host plant defense

A B S T R A C T Cotton (Gossypium hirsutum) is an important fiber crop worldwide. Insect attack causes cotton yield and quality losses. However, little is known about the mechanism of cotton response to insect attack. We simulated insect feeding by applying insect oral secretions(OS)to wounds,and combined transcriptome and metabolome analysis to investigate how OS from two major pest species (Helicoverpa armigera and Spodoptera litura) affect cotton defense responses. We found that respectively 12,668 and 13,379 genes were differentially expressed in comparison with wounding alone. On addition of OS, the jasmonic acid signaling pathway was rapidly and strongly induced, whereas genes involved in salicylic acid biosynthesis were downregulated. On constructing a coexpression gene network, we identified a hub gene encoding a leucine-rich repeat receptor kinase that may play an important role in early signal recognition and transduction.OS from the two insect species altered the abundance of flavonoid-related compounds in different patterns. Gossypol remained in lower concentration after OS application than after wounding alone,suggesting a suppressive effect of OS on cotton defense response. This study illustrated transcriptional and metabolic changes of cotton in responding to OS from two chewing insect species, identified potential key response genes, and revealed evidence for OS inhibition of wounding-induced cotton defense response.

1.Introduction

Plants defend themselves against herbivories using various strategies.Evidence has accumulated over the past few years to indicate that hormones,via the jasmonic acid(JA),salicylic acid (SA), and ethylene (ET) signaling pathways, play central roles in regulation of plant defense against pathogen attack and herbivory [1–4]. Transmembrane receptor-like protein kinases (RLKs), among the most important sensory protein groups for plant perception of environmental cues such as pathogen-associated molecular patterns (PAMPs) or herbivore-associated molecular patterns (HAMPs), regulate plant defense response to pathogens and herbivores[5–7].

Plants also produce structurally and functionally diverse secondary metabolites, of which toxic or anti-digestive compounds may inhibit feeding [8–11]. These include benzoxazinoids and phenolamides in Poaceae or glucosinolates in Brassicaceae [12,13]. In cotton, compounds including rutin and gossypol showed anti-insect activity[14].

In turn,herbivores have developed strategies to counter plant resistance. The most straightforward strategy is elevated tolerance to the defensive compounds synthesized by host plants.For example, a CYP6AE gene cluster in the Helicoverpa armigera genome contributes to detoxification of phytochemicals and insecticides [15]. Some herbivores release specific effector molecules or insect-associated bacteria that suppress defense response in host plants[16–20].Insect endosymbiotic bacteria may influence the interaction by expanding the food range or changing the phenotype of plants to benefit their insect attackers[21,22].

Unlike abiotic stresses,the feeding behavior of herbivores,such as feeding time and feeding site,is difficult to control.A high degree of treatment standardization is needed for the investigation of mechanisms. The application of oral secretions(OS)to mechanical wounds[23,24]has been found[25]to mimic herbivore feeding effectively.

Cotton (Gossypium hirsutum) is an important worldwide source of natural fiber and the third largest source of edible oilseed tonnage [26,27]. However, cotton production is challenged worldwide by diverse arthropod pests, occurring at economically damaging levels every year [28]. The herbivorous insects cotton bollworm (H. armigera) and cotton leafworm (Spodoptera litura) are the two most destructive herbivores feeding on cotton. They feed on the cotton leaf,bud and boll, resulting in damage to cotton growth and production. Although the interactions between insects and host plants have been studied intensively [13,29–31], little is known about the defense response of cotton against chewing insects at the transcription and metabolic levels.

In this study, we investigated the interaction between these two chewing insects and cotton plants by mimicking insect feeding. We performed transcriptomic and untargeted metabolomic analyses to identify differentially responding genes and metabolites in cotton leaves subjected to wounding and OS treatments.

2. Materials and methods

2.1. Plant treatment and sample preparation

The cotton cultivar Jin 668 was sown in nutrient-rich soil in plastic pots with controlled temperature and daylight hours(28–30°C,16 h).Second-instar larvae of H.armigera and S.litura were purchased from Henan Jiyuan Baiyun Industry Co.,Ltd.,Henan, China. Insects were reared in 9-cm Petri dishes containing fresh leaves of Jin 668.

OSs were collected from larvae(third to fifth instar)and kept on ice during collection and then stored at ?70 °C until use.Plants used for the experiment were 30 days old and the third leaf from the top of each plant was selected for treatment.Wound+water(W+W)and wound+OS(W+OS)treatments were performed by application of 10 μL of water or OS to wounds generated by rolling a sewing pattern wheel along the midvein to generate two rows of wounds on each side of two fully expanded leaves per plant (Fig. S1-a–f). Treatments were all performed at the indicated times prior to harvest. Treated leaves from three individual cotton plants were harvested as one experimental replicate,flash-frozen in liquid nitrogen,and stored at ?80°C.For each time point,respectively three and six replicates were collected for RNA and metabolite extraction.

2.2. RNA isolation and library construction

To assess the global transcriptome profile of cotton plants regulated by W + OS and W + W treatment, we performed deep RNA sequencing of the leaf samples treated for 0,10,30,60, and 240 min, with untreated plants indicated as 0 min.Total RNA was extracted using a polysaccharide and polyphenolic-rich RNApre Pure Plant Plus Kit (Tiangen Biotech,Beijing,China).RNA concentrations were measured with a Qubit 2.0 fluorometer (Invitrogen, CA, USA) and RNA integrity was assessed using Agilent Bioanalyzer 2100 system(Agilent,CA,USA).RNA samples of 1.5 μg were used for library construction. RNA sequencing was performed on an Illumina HiSeq 4000 platform(Illumina, CA,USA).

2.3. Analysis of RNA-sequencing data

Clean reads were obtained by filtering out low-quality reads.Mapping, quantification of expression, and differential expression analysis were performed with TopHat,Cufflinks and Cuffdiff [32,33]. Genes with P (adjusted) < 0.05 and | log2(fold change) | ≥1 were defined as differentially expressed genes(DEGs). The genome of G. hirsutum cv. TM-1 was used as the reference genome for reads mapping [34]. Gene expression profile data for Pearson correlation coefficient matrix and principal component analysis (PCA) plots were preprocessed by Pareto scaling.Gene Ontology terms were assigned to DEGs using WEGO 2.0 [35]. Metabolic pathways were predicted by Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis using KOBAS 3.0[36].

2.4. qRT-PCR analyses

To further confirm the expression level of genes, qRT-PCR analyses was performed. First-strand cDNA was generated from 2 mg of total RNA using SuperScript III reverse transcriptase(Invitrogen,CA,USA).qRT-PCR was performed in 15 mL reactions using the ABI 7500 Real-time PCR System (ABI, CA, USA), with three technical replicates and three independent biological replicates. GhUB7 was used as an internal control. The primers used in our study are searched from qPrimerDB[57].

2.5. Measurement of phytohormones

Cotton leaves (100 mg fresh weight) were homogenized with 80%cold methanol(v/v)with internal standard(10 ng mL?1of 9,10-dihydro-JA)and shaken at 4°C overnight in the dark.The extract was filtered through a 0.22-mm nylon membrane and stored at ?70 °C prior to measurement. Quantification of the phytohormones was performed on an Applied Biosystems 4000 QTrap high-performance liquid chromatography-mass spectrometry system(ABI,CA,USA).

2.6. Untargeted metabolomics data acquisition

Secondary metabolites are considered to play important roles in plant adaption to environment stress. To find metabolites related to cotton defense response induced by herbivory, an untargeted metabolomics experiment was performed(Fig.S4-a). We analyzed the content of metabolites extracted from cotton leaves which were induced by W + W and W + OS treatment after 1 h,4 h,10 h,and manually annotated some of them(Fig.S4-b).In the untargeted metabolomics experiment,a Waters Acquity UPLC system (Waters, MA, USA) equipped with a BEH C18 column(1.7 μm,2.1 mm×100 mm)was used to separate plant extracts using a gradient program with a flow rate of 400 μL min?1. The gradient program was as follows:0–1 min,99%A,1%B;1–2 min gradient phase to reach 80%A,20%B;2–6 min,60%A,40%B;6–7 min,45%A,55%B;7–15 min,5% A, 95% B; 15–18 min, 5% A, 95% B (solvent A: 99.0% water,1.0% methanol, 0.1% formic acid; solvent B: 5.0% water, 95%methanol, 0.1% formic acid. All percentages are v/v). The temperatures of the column and sample chamber were held at 45°C and 4°C,respectively.A Waters Xevo G2xs Q-TOF MS fitted with an ESI source was used to acquire MS and MSndata. Data were collected in both positive and negative ion mode by the MSE method,using respectively 6 V and 15–30 V collision energies for the low- and high-energy channels.Calibration was performed prior to sample analysis via infusion of a sodium formate solution with mass accuracy within 1 μmol L?1. The capillary voltage was set to 2500 V for positive ion mode,the source temperature was held at 100°C,and the desolvation temperature was set to 400 °C with a nitrogen desolvation gas flow rate of 800 L h?1.

2.7. Metabolomic data processing

Raw data files from the UPLC-MSE platform were converted to CDF format with Databridge software (Waters, MA, USA).Those converted files were further processed with XCMS [37]and CAMERA [38] (all based on R 3.3.0 software), including peak detection, retention time correction, peak filling, and annotation of isotopes and adducts.The “centWave” method was selected to perform peak detection, using parameters ppm = 20, snthresh = 10, peakwidth = c (2, 10). For retention time correction,the parameters were set as follows:minfrac=0.7, bw1 = 3, bw2 = 2, mzwid = 0.015. Isotope annotation and in-source pseudospectra were reconstructed with CAMERA.The LOESS method-based quality control (QC) samples [39]were used to perform correction for each metabolite signal(feature) detected in our experiment and only features with relative standard deviation (RSD) < 30% in QC samples were retained for further analysis. The dataset was submitted to MetaboAnalyst 3.0 [40] for statistics analysis. Features with> 50% missing values were removed, and the remaining missing values were replaced by a small value (half the minimum positive value in the original data). To improve the quality of data analysis, the data set was filtered by interquartile range (IQR) removing uninformative variables that showed near-constant values throughout the experiment.Features with P<0.05(t-test,false discovery rate(FDR)adjusted),and fold change<0.67 or>1.50 were identified as differential features.

2.8. Metabolite identification

The molecular formulas of features were predicted by MassLynx software(Waters,MA,USA)based on accurate m/z and isotopic pattern in positive model (window < 10 μmol L?1). Mass fragments of putative molecules were compared with those of candidate compounds in online mass spectral libraries(METLIN,Massbank,HMDB).Some of them were identified by comparison with spectra of purchased standards.

2.9. Data availability

The transcriptome data set generated during the current study has been submitted to NCBI under accession number PRJNA522889. The scripts used for metabolomics data processing can be downloaded from https://github.com/huans/Working-scripts.

3. Results

3.1. RNA sequencing of cotton treated with mechanical wounding and oral secretions (OS) from two chewing insects

In total, 39 cDNA libraries were constructed, yielding 288 Gb of sequence with ~25 million reads per library, of which ~1% of paired-end reads were filtered and trimmed (Table 1). The correlation coefficients between pairs of biological replicates under the same treatment were all greater than 0.95,indicating satisfactory quality of the RNA sequencing (RNA-Seq) data (Fig.S2).

3.2.Transcriptome profiling of cotton leaves after H.armigera OS application

A total of 12,668 DEGs were detected in W+W and W+OS_H treated samples at least one time point (Table S1). In this analysis, the 10-, 30-, and 60-min time points of the W + OS samples clustered farthest from W+W,whereas samples in 240-min clustered close together, indicating that gene expression level changed rapidly after OS application and decayed in a short time (Fig. 1-a). This inference is supported by the numbers of DEGs at each time point(Fig.1-b).Plants responded to W+OS_H treatment by differentially expressing 5494 genes at 10 min,increasing to 8153 at 30 min,reaching the maximum of 8495 at 60 min and decreasing to 340 at 240 min. In total 1164 upregulated and 1193 downregulated genes were observed at the first three time points, whereas only 13 and 58 DEGs were observed at all time points(Fig.1-c,d).

DEGs with upregulation were assigned to 17 KEGG pathways(P < 0.05), including plant-pathogen interaction, alpha-linolenic acid metabolism,and plant hormone signal transduction.DEGs with downregulation were assigned to photosynthesis, carbon metabolism and 39 other pathways(Fig.1-e,Table S2).

3.3.Transcriptome profiling of cotton leaves after S.litura OS application

Consistent with OS_H,plants treated by OS_S exhibited rapid response and decayed in a short time(Fig.2-a,b).In comparisonwith W + W, a total of 13,379 DEGs were found for W + OS_S,consisting of respectively 5130,7747,8894,and 330 genes at the 10-, 30-, 60-, and 240-min time points (Fig. 2-b, Table S3). In a pattern similar to that of W+OS_H treatment,W+OS_S-treated samples also showed a larger proportion of DEGS with downregulation than with upregulation. More than 6000 genes showed significant responses to OS in a short time.Most of the DEGs were observed for the first three time points,with 3253 consistently upor downregulated at all time points except for 240 min(Fig.2-c,d),a response consistent with that to W + OS_H. KEGG pathway analysis assigned all up-or downregulated DEGs to respectively 23 and 39 pathways(P<0.05)(Table S4).Upregulated genes were enriched in alpha-linolenic acid metabolism and plant-pathogen interaction and downregulated genes were enriched in photosynthesis(Fig.2-e).

Table 1–Description of RNA sequence from water-or OS-treated Jin 668 leaves.

Comparison of the DEGs from the two OS treatments showed that respectively 4471 and 5433 genes were commonly up- or downregulated after treatment with W + OS_H and W+OS_S(Fig.S3-a,b).The GO analysis showed that most of the terms to which DEGs were assigned were common to both OS treatments (Fig. S3-c). DEGs were significantly enriched in catalytic activity (GO:0003824), binding(GO:0005488), metabolic process (GO:0008152), and cellular process (GO:0009987). Terms of genes involved in defense response to bacterium and fungus were also enriched significantly (GO:0050832,GO:0042742, GO:0009615).

3.4. OS induced early signaling pathways

Most major plant responses to herbivores or pathogen attack are regulated by JA,SA,and ET.KEGG pathway analysis of DEGs(Figs.1-e,2-e)suggested that cotton defense signaling pathways are significantly influenced by OS. JA-, SA-, and ET-responsive genes include lipoxigenase (LOX), allene oxide synthase (AOS),allene oxide cyclase (AOC), oxophytodienoate-reductase 3(OPR3), 1-aminocyclopropane-1-carboxylate synthase 1 (ACS1),1-aminocyclopropane-1-carboxylic oxidase (ACO). LOX2 was upregulated after the onset of OS application, and AOC, AOS,and OPR3 were significantly upregulated after 30 min and reached their highest levels at 60 min(Fig.3-a),as confirmed by qRT-PCR(Fig.3-b).

Fig.1–Overview of cotton transcriptome response to oral secretions from H.armigera.(a)PCA plot of genes identified by RNASeq of inbred line Jin 668 treated with W+W and W+OS for the indicated times.(b)Number of DEGs at each time point.P<0.05(FDR-adjusted), and fold change< 0.5 or >2.0. (c,d)Venn diagram showing the distribution and overlap of up- and downregulated DEGs across all time points.(e)The top seven pathways enriched in significantly up-(left) or down-(right)regulated genes(P <0.05). Circle size represents the number of enriched genes.The X axis displays the enrichment factor.

Expression of vegetative storage protein 2 (VSP2), an important JA-responsive gene, also significantly increased at 60 min(Fig.3-a).This result indicated that OS from both insect species induced the expression of JA-related genes.JA and JAIle rapidly increased in both W + W and W + OS relative to control, and the levels induced by W + OS_H and W + OS_S were more than one-fold greater than induced by W+W(Fig.3-c), in agreement with the changes at the transcriptional level.

In contrast,expression of genes involved in SA biosynthesis, phenylalanine ammonia-lyase (PAL) and isochorismate synthase (ICS), decreased significantly in W + OS in comparison with W + W (Fig. 3-a). The expression of pathogenesisrelated gene (PR1, PR5), a SA-responsive gene, showed no significant change (data not shown). Ten minutes after wounding,the transcripts of PAL and ICS rapidly accumulated and maintained a high level across all time points in W + W relative to control (non-treatment) (Fig. 3-a), suggesting that they could be induced by wounding and suppressed by OS.Although the ET-biosynthesis related genes ACS1 and ACO were induced by OS, expression of the ET-responsive genes ethylene-insensitive 3 and ethylene response factor 1 (EIN3,ERF1) decreased significantly at 30 min and 60 min (Fig. 3-a),indicating an underlying mechanism in ethylene signal transduction during plant-insect interaction.

3.5. Co-expression network helps to identify defense response related hub gene

To identify key genes(hub genes)that respond to OS,a weighted correlation network analysis of the DEGs was performed(Fig.4-a). The expression matrix of transcripts from 39 samples was used for this analysis,with samples treated with W+W at 10,30,and 60 min excluded as outliers. The analysis identified 10 significant co-expression modules (Fig. 4-b). The 19 genes that could not be clustered into any module were assigned to a gray module.KEGG enrichment analysis was performed for genes in each module (Table 2). According to the results of enrichment analysis,the blue module,which showed significant association with plant-pathogen interaction, was selected for further analysis(Fig.4-c).

Fig.2–Overview of cotton transcriptome response to oral secretions from S.litura.(a)PCA plot of genes identified by RNA-Seq of inbred line Jin 668 treated with W+ W and W+OS for the indicated times.(b) Number of DEGs at each time point.P <0.05(FDR-adjusted), and fold change< 0.5 or >2.0. (c,d)Venn diagram showing the distribution and overlap of up- and downregulated DEGs across all time points.(e)The top seven pathways enriched in significantly up-(left) or down-(right)regulated genes(P <0.05).Circle size represents the number of enriched genes.The X axis displays the enrichment factor.

To investigate how genes in this module responded to OS,we examined the expression profiles of genes involved in plant-pathogen interaction (Fig. 4-d). We found that they showed similar expression patterns in W + OS_H and W + OS_S, with most of them upregulated after treatment with OS in comparison with W + W at the same stage.Examples are kinase heat shock protein(HST)and respiratory burst oxidase homolog a (RBOHA) and WRKY22, which are involved in defense response. The abscisic acid (ABA) signaling genes calcium-dependent protein 4 (CPK4, Gh_A03G1505,Gh_A13G0566, Gh_D02G1973) and calcium-dependent protein 32 (CPK32, Gh_D09G1163), which regulate ABA-responsive gene expression by phosphorylating ABF4, were significantly induced at 10,30,and 60 min after OS application (Fig.4-d).

The search for hub genes in this module was based on high connectivity with other genes. Gh_D10G0074, which has high amino acid sequence identity(55%)with pep1 receptor 1(PEPR 1), encoding a leucine-rich repeat receptor kinase in Arabidopsis thaliana, was identified as the hub. The transcript abundance of Gh_D10G0074 was induced rapidly at 10 min after OS application and at least 3-fold upregulated relative to W +W at 30 min(Fig.4-e).

3.6. Overview of metabolic profiling of cotton response to oral secretions

After data preprocessing and filtering, a total of 11,155 features were remained, and 9809 with RSD < 0.3 in QCs were selected for further analysis (Fig. S4-c). We used all samples (including QCs and subject samples) for principal component analysis, finding that all of the QC samples were clustered well(Fig.S4-d).

Like the transcriptome profiles, metabolomic profiles of cotton leaves treated with W + OS showed sharp differences compared with W + W. PCA with all filtered metabolite features showed a clear separation between W + OS and W + W treatment especially at 1 h and 4 h (Fig. 5-a, b),suggesting that there were large metabolic differences between plant response to mechanical wounding and OS. In addition, the samples at 1 h with W + OS_H treatment clustered farthest from W + W at the same stage, whereas 10-h W + OS_H and W + W sample remained relatively close and were closer to control (untreated plants) at other stages(Fig.5-a),indicating that metabolic differences disappeared as time increased. However, a different pattern was observed in W + OS_H, which showed clear differences even at the last time points(Fig.5-b).

Fig.3– Expression of genes and metabolites involved in hormones.(a)Expression profiles of JA-(indicated by red line), ET-(blue line)and SA-(green line)related biosynthesis or responsive genes.(b)Validation of the gene expression by qRT-RCR.(c)Abundance of JA and JA-Ile after different treatments.Values are mean± SE(n =5). *,P <0.05 by t-test relative to W+ W.

Table 2–KEGG enrichment pathway for co-expression modules.

Respectively 297, 330, and 373 features were significantly upregulated and 1305,169,and 171 were down-regulated after 1 h, 4 h, and 10 h (P (FDR adjusted) < 0.05, with fold change >1.50 or<0.67)respectively in W+OS_H relative to W+W(Fig.5-c).At the initiation of W+OS_S treatment,there were many more features with up- than downregulation, whereas the reverse was seen at 4 h, with a larger number of downregulated features observed at 10 after W+OS_S treatment(Fig.5-d).There was limited overlap across the three time points(Fig.5-e, f), indicating the dynamic changes in cotton plant response to OS.

3.7. OS-induced changes in flavonoid-related metabolism in cotton plants

Flavonoids in plants are involved in a multitude of functions including as pigments and antioxidants, with some of them participating in plant defense against biotic and abiotic stresses[41,42]. To investigate the effects at the metabolomic level in cotton response to OS, we focused on several compounds annotated as flavonoids or flavonoid-related(Fig.6-a).

Fig.5– Overview of cotton metabolic response to W+W and W+ OS.PCA plot of metabolite features of inbred line Jin 668 treated with W+W,W+OS_H(a)or W+OS_S(b)for the indicated times.Number of significantly regulated metabolite features in W+ OS_H(c)and W+OS_S(d)respectively at each time point.Venn diagram showing the distribution and overlap of significantly changed metabolite features in W+ OS_H(e)and W+ OS_S(f). P <0.05 (FDR adjusted), and fold change<0.67 or> 1.50.

Fig.6– Effects of OS and wounding on flavonoid-related metabolites.(a) Flavonoid-related biosynthesis pathway. (b–i)Metabolites that were significantly altered in content by OS.Values are mean± SE(n = 3).

In our study, we found that chlorogenic acid, syringin, and procyanidin B5 could be induced after treating by W+OS_H for at least one time point (Fig. 6-b–d). Interestingly, the content of syringin and procyanidin B5 in W + OS were downregulated at 10 h in comparison with W + W (Fig. 6-b–d), suggesting an inhibitory effect of OS on plant defense response. Isoquercitrin and rutin, which have been shown [14] to be important in protecting cotton from herbivore feeding, did not show significant differences between W + W and W + OS. For other compounds, quercetin 3-O-malonylglucoside and quercetin 3-xylosyl-(1->2)-alpha-L-arabinofuranoside, derived from isoquercitrin and quercetin respectively, showed significant increases at 10 h after W + OS_H, whereas no significant difference was observed in W+OS_S(Fig.6-g,i).The contents of quercetin 3-sambubioside and quercetin 3-(2G-xylosylrutinoside)were also significantly induced at 10 h, and kaempferol 3-Oglucoside was observed to increase at the last two stages in W + OS_H (Fig. 6-f, h). Induction effects for these three compounds were found in W + OS_S at 10 h, but they were not statistically significant except for kaempferol 3-O-glucoside.These results indicated that OS from the two insect species affected the accumulation of flavonoids with different patterns,and it appeared that cotton plants responded more strongly to OS from H.armigera.

3.8. Gossypol biosynthesis pathway was influenced by OS

As a representative sesquiterpene derivative compound in cotton,gossypol is considered[14]to play an important role in cotton against insect infestation. To investigate changes in gossypol,we analyzed the expression of genes involved in the gossypol biosynthesis pathway,including(+)-δ-cadinene synthase (CDN), cytochrome p450 monooxygenase, (CYP706B1,CYP82D113, CYP71BE79), alcohol dehydrogenase (DH1), and 2-oxoglutarate/Fe(ii)-dependent dioxygenase(2-ODD-1)(Fig.7).

Almost all of the genes involved in this pathway showed higher expression levels at 4 h after W+W treatment;once OS was applied, the abundance of transcripts was significantly reduced (Fig. 7). These results suggest that the OS from the two insect species suppressed the gossypol biosynthesis by downregulating the expression of genes in this pathway. To test this hypothesis, we measured the content of gossypol after treatment with W+W and W+OS at several time points.The content of gossypol reached its highest level at 4 h after wounding, and remained in lower abundance after W + OS treatments throughout the experiment, consistent with the pattern at the transcriptional level (Fig. 7). These results suggested that some specific components in OS play important roles in plant-insect interaction by influencing plant defense response.

Fig.7–Expression of genes in the gossypol biosynthesis pathway and the abundance of gossypol after treatment with W+W and W+ OS.Values are mean±SE (n=5).**,P< 0.01 by t-test relative to W+ W.CNDC,(+)-δ-cadinene synthase;CYP,cytochrome p450 monooxygenase;DH1,alcohol dehydrogenase;2-ODD-1,2-oxoglutarate/Fe(ii)-dependent dioxygenase.

4.Discussion

Oral secretions from herbivores influence plant defense response [43,44]. Although only a very small amount of OS(10 μL)was applied to artificial wounds,significant differences in transcriptional and metabolic responses of cotton to W + OS could be observed in a very short time. However,transcriptional differences almost vanished at 240 min (Figs.1-b, 2-b), whereas in a similar study in maize [23], more than 2000 genes still showed significant regulation 6 h after treatment with Mythimna separata OS, possibly owing to the species specificity of host plants and insects.

In our study, significant overlap was observed at the transcriptome level after treatment with W + OS_H and W + OS_S (Fig. S3), indicating the similar response patterns of cotton to H. armigera and S. litura attack. KEGG pathway enrichment analysis showed that genes involved in photosynthetic pathways were significantly downregulated after treatment with OS (Figs. 1-e, 2-e). This finding is consistent with the concept that pathogen infection can inhibit photosynthesis, including photosystem II (PSII) activity, CO2fixation,and globally downregulate photosynthetic genes[45,46].In a recent study [47], effector-triggered immunity (ETI)inhibited photosynthesis by activating mitogen-activated protein kinases 3/6 (MPK3/6). In our study, MAPK3 was induced rapidly after OS treatment (Fig. S5). These results suggest that in our study OS may also have induced cotton plant ETI as do pathogens.

As a major hormone in plant signaling, JA plays an important role in plant defense against chewing insects and insect OS[48,49].Mutant plants with deficiency in JA signaling are more susceptible to insects than wild-type plants [50,51].In our study,the contents of JA or JA-ILE increased after cotton plants were treated with W + OS, in comparison with plants treated with W+W at all time points(Fig.3-c).With respect to another important signaling pathway, SA-biosynthesis related genes were significantly downregulated (Fig. 3-a). The interaction between ABA and ET signaling is mutually antagonistic in vegetative tissues [52]. This relationship may explain why ET-biosynthesis genes were upregulated but ETregulated genes such as EIN3 and ERF1, were significantly downregulated, as CDPK4 and CDPK32, which positively regulate ABA signaling, were strongly induced (Fig. 4-d).These results indicate that the JA and ABA signaling pathways were dominant in cotton response to OS.

Co-expression network analysis revealed a hub gene(GhD10G0074) in the module specifically enriched in the plant-pathogen interaction pathway. It has high amino acid sequence identity (55%) with PEPR 1 encoding a leucine-rich repeat receptor kinase in Arabidopsis thaliana. In previous studies, PEPR 1 has been reported as a receptor for endogenous defense signal, Pep1, a 23-aa peptide that activates transcription of a defensive gene (PDF1.2) and activates the synthesis of H2O2to amplify innate immune response to pathogen attacks[53–56].In our study,the expression of PEPR 1 was significantly induced at the outset of OS treatment(Fig.4-e), indicating the rapid recognition of cotton to signals.These results also demonstrate that early signal recognition and transduction are conserved in higher plants.

Plants produce a large number of functionally diverse secondary metabolites, some of which play important roles in plant response to biotic stress [11]. In our study, OS from two chewing species affected the accumulation of flavonoids in different patterns, but OS from S. litura was weaker than OS from H. armigera. OS from H. armigera induced the accumulation of some anti-insect substances such as syringin and procyanidin B5, but OS from S. litura did not, and even slightly inhibited this process. As one of the most important compounds toxic to herbivores, gossypol was downregulated after treatment with W + OS at 4 h in comparison with W + W (Fig.7). These results remind us that as a specific effector or microbe in OS may play a critical role in plant-insect interaction by influencing plant defense response.

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2019.12.007.

Author contributions

SX and XL designed and coordinated the study. HS, HL, YW,ZX, SL, BL, XD, JL, LS, QW carried out all the experiments. HS analyzed transcriptome and metabolomics data and wrote the draft manuscript. All authors approved the final manuscript.

Declaration of competing interest

All the authors have declared no conflict of interest.

Acknowledgments

This work was supported by the National Key Research and Development Program of China (2016YFD0100203-9), National R&D Project of Transgenic Crops(2016ZX08010001-006),National Natural Science Foundation of China (31371673). We thank Dr.Jianqiang Wu at Yunnan Key Laboratory for Wild Plant Resources, Kunming Institute of Botany, Chinese Academy of Sciences for providing help in manuscript writing.