Comparative transcriptome analysis between thrips-resistant and thrips-susceptible alfalfa (Medicago sativa L.)

Background: Thrips (Thysanoptera: Thripidae) are major insect pests on alfalfa and result in decreased plant nutrients and growth, low yields and even plant death. In our previous studies, an alfalfa variety (Caoyuan No.4) with high thrips resistance was bred through consecutive eld recurrent selection. In order to better understand the genetic and molecular mechanisms of thrips resistance in Caoyuan No.4, RNA-Sequencing was employed using the thrips-resistant alfalfa accession (Caoyuan No.4) and a thrips-susceptible alfalfa accession (Caoyuan No.2), each with and without thrips infestation. Results: There were 851 genes constitutively upregulated and 434 genes downregulated in Caoyuan No.4 compared to Caoyuan No.2 without thrips infestation. The upregulated genes were mainly involved in primary metabolism such as energy metabolism and carbohydrate metabolism, lipid metabolism and certain secondary metabolites, while the downregulated genes were mainly related to plant-pathogen interaction. In addition, very few DEGs (only 13) were detected in Caoyuan No.4 after thrips stress, but a total of 3326 contigs DEGs were detected in Caoyuan No.2 after thrips stress. The upregulated genes in Caoyuan No.2 after stress were mainly involved in isoavonoid biosynthesis, proteasome, amino sugar and nucleotide sugar metabolism, avonoid biosynthesis as well as plant-pathogen interaction. Moreover, 117 genes that were shared in both the S_CK vs S_T group and S_CK vs R_CK group were divided into 6 clusters, which are mainly involved in secondary metabolism, fatty acid metabolism, amino acid metabolism, rust resistance kinase, WRKY transcription factor and nodule lectin. Conclusion: Both constitutive defensive genes and potential induced defensive genes were detected in the defense of Caoyuan No.4. That two distinct kinds of defensive genes — constitutive defensive genes and induced defensive genes — can be simultaneously activated and thus potentially enhance plant protection against insects attacks is a signi ﬁ cant nding for plant resistance breeders.

associated with host plant resistance to thrips such as leaf hair density, leaf hairiness, leaf hardness, leaf wax, glandular hairs, and trichomes [8][9][10][11]. However, biochemical-based defence is considered more effective as it directly affects insect growth and development [12]. It has been reported that many secondary metabolites and plant hormones are associated with thrips resistance, such as protease inhibitors, phenols, tannins, salicylic acid (SA) and jasmonic acid (JA) [13][14][15][16] Since there is an increasing focus on improving crop production through safe and sustainable means by reducing the reliance on pesticides [17], the use of resistant cultivars is currently considered the most effective and environmentally sustainable strategy to control insects. Alfalfa cultivars are heterogeneous populations, and recurrent selection is highly heritable and it is relatively easy to accumulate excellent genes through a large population. Generalist and specialist insect herbivores show an evolution of some candidate genes responsible for their adaptation to host plants [18]. Eco-genomic tools have been used to study the genetic basis of plant defensive traits in many plants [19]. However, little is known about the molecular and genetic mechanisms involved in alfalfa thrips resistance. Next generation sequencing and more speci cally RNA-Sequencing (RNA-Seq) has become a popular and comprehensively informative approach to monitor wide transcriptional changes during insect-plant co-evolution. Recently, a comparative transcriptomic analysis was employed to assemble the expressed transcriptomes of alfalfa, which mainly focused on the induced defence genes related to both resistant and susceptible alfalfa lines after thrips infestation [3].
In our previous study, an alfalfa variety (Caoyuan No.4 ) with high thrips resistance (hazard point coe cient: 0.26, and pest index: 0.33) was bred through nearly 30 years of consecutive eld recurrent selection, and the mechanisms of thrips resistance in Caoyuan No.4 were preliminarily studied anatomically and biochemically [20,21]. However, the genetic and molecular mechanisms underlying Caoyuan No.4 resistance to thrips are not well understood. Based on previous anatomical and biochemical studies [20], we hypothesized that the molecular mechanisms supporting constitutive defence may be more important for thrips resistance in Caoyuan No. 4 [21].

Alfalfa treatement
Alfalfa plants were treated as described by Xiong et al. [3] with some modi cations. When the seedlings reached budding stage (about 60 days), both cultivars were randomly and equally divided into two groups: (1) 30 alfalfa thrips (Odontothrips loti) per plant were placed onto the leaves and covered by a cage with a 90-mesh nylon cloth, resulting in the R_T (No.4) and S_T (No.2) treatment groups; and (2) plants were not treated with thrips and were maintained under the same conditions, which led to plants in the R_CK (No.4) and S_CK (No.2) treatment groups. After 3 weeks, the top 3-4 leaves were cleaned of any thrips and harvested from each treatment. All samples were immediately frozen in liquid nitrogen and stored at -80 o C.
RNA extraction, cDNA library construction and RNA-sequencing Three biological replicates were used for all RNA-Seq experiments from each cultivar and thrips treatment. RNA extraction, cDNA library construction and RNA-sequencing were carried out as described by [22,23]. Brie y, total RNA was extracted from the leaves using Trizol reagent (Invitrogen, Carlsbad, CA) along with DNase treatment, according to the manufacturer's instructions (QIAGEN, Germany). The quality and quantity of total RNA were assessed using NanoDrop 2000 analysis and gel electrophoresis. As described by [24], the cDNA libraries were prepared using a TruseqTM RNA sample prep Kit (Illumina), and RNA-sequencing was performed on an Illumina Hiseq 4000 (Version 2 × 150 bp) at Shanghai Majorbio Bio-pharm Biotechnology Co., Ltd. (Shanghai, China). The raw sequence reads were deposited in the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/Traces/sra, accession number PRJNA622603).

De Novo Assembly
The clean data were obtained by removing the adapter and primer sequences using SeqPrep software, and fragments of less than 20 bp in length were excluded from further analyses using software Trinity (https://github.com/trinityrnaseq/trinityrnaseq) in the absence of a reference genom [25]. The sequence assembly quality was evaluated using the number of sequences and bases, GC percentage, distribution of unigene lengths, average coverage, and N50 statistics [26].

Annotation and Classi cation
The assembled transcriptome sequences were searched against six databases (NR, Swiss-Prot, Pfam, COG, GO and KEGG databases) to obtain annotation information in each database. Speci cally, to obtain the similarity to other species and the functional information of homologous sequences, the sequences were searched against the NCBI non-redundant database (NR, ftp://ftp.ncbi.nlm.nih.gov/blast/db/), Swiss-Prot database (http://web.expasy.org/docs/swiss-prot_guideline.html) and Pfam (http://pfam.xfam.org/). The gene function terms were obtained through the Gene Ontology database (GO, http://www.geneontology.org). Functional classi cation was performed using the Clusters of Orthologous Groups of proteins database (COG, http://www.ncbi.nlm.nih.gov/COG/), and pathway annotation was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/). A P value < = 0.05 was regarded as the threshold for signi cance [26]. In addition, quantitative analysis of gene and transcript expression levels were obtained through RSEM (http://deweylab.github.io/RSEM/), and principal component analysis (PCA) was performed to obtain the relationships among and variability between samples .

Differentially Expressed Gene Analysis And Functional Enrichment
Differential expression gene analysis of the samples was performed using the DEseq 2 software. A padjust value < 0.01 & |log2FC | >=1 was set as the threshold. To identify shared and unique genes/transcripts across gene sets, venn analysis was performed. Gene set enrichment analysis of the differentially expressed genes (DEGs) was then performed for the KEGG annotations to determine overrepresented functional pathways (with a p-value < 0.05) at each comparison level (for different genotype and treatment) [27].

RT-qPCR Validation
Total RNA was extracted from leaves with different treatments and rst strand cDNA was synthesized with a PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Dalian, China) according to the manufacturer's instructions. The qRT-PCR was performed with three biological replicates and three technical replicates as previously described [28]. The genes and primer sequences used in qRT-PCR assay are shown in Table 1. Table 1 The primers for qRT-PCR in this study TAGCCATCCCGAAATTCGCT  TCGTTGGTGCCATCCAGTTT  TTGCACCCAACTCATAGCCC  ACTATGAAAACAAAGGGTGGG  T  GATACAAAGTCAAAGCACTATC  ACA  TAAAACCGTGCTGGAACCCT  GAAGCACAAGAGTGGCTCC  ACAGAATCGTCGACAAGGCA  TCTGGAACCAGCCATTGACC  TATGGGTGGAACGAGACCCA  TGGTCGAATTCCAGGCCTTT  TTCCCGAGCTTGTAGTGTGC  TGCCACTAAACGAGGTCCAG  TTGCCGAAGGTAAGAACACG  AACGAAGGCAAGTCCCCAAT  ATCACCAGCACCATCAGGAC  TGTGAGTTGATTCGGTCGGA (Table S1). A transcriptome datebase containing 576 988 transcripts of average length 607.22 bp was obtained using Trinity software, with an N50 length of 883 bp and an E90N50 length of 1303 bp. Moreover, the sequence length of the assembled transcripts in 61% of sequences was < 500 bp in length with a minimum of 201 bp and a maximum of 15,749 bp ( Table 2).  3 . Therefore, this sample was excluded from subsequent DEGs analysis (Fig. 2).
The DEGs among the four sampling groups were analyzed in the RNA-seq datasets. As shown in Fig. 3, a total of 5,655 signi cant DEGs were detected. Constitutive differential expression was shown in the control group (S_CK vs R_CK), with 851 contigs constitutively upregulated and 434 contigs downregulated in R_CK. After treatment with thrips, a total of 2,865 upregulated and 461 downregulated contigs were detected in S_T group (S_CK vs S_T), while only 11 and 2 upregulated and downregulated contigs were found in Caoyuan No.4 alfalfa (R_CK vs R_T). In addition, 2260 and 1064 upregulated and downregulated contigs were observed between the two groups treated with thrips (R_T vs S_T) (p-adjust value < 0.01 & |log 2 FC | >=1) (Fig. 3).
The result showed that 117 genes were shared in both the S_CK vs S_T group and the S_CK vs R_CK group, while the other 2087 DEGs were only observed in the S_CK vs R_CK group (Fig. 5a). We hypothesize that the shared genes in the S_CK vs S_T group might be the induced defensive genes, while the other 2087 DEGs might be the constitutive defensive genes to thrips in Caoyuan No.4. A heat map generated using the 117 shared DEGs was divided into 6 clusters (Fig. 5b, Table S4). For cluster 1, all the 27 DEGs in both S_CK and S_T were upregulated, while the genes of R_CK and R_T were downregulated. These genes were mainly enriched in stilbenoid, diarylheptanoid and gingerol biosynthesis. Cluster 2 contained 22 DEGs upregulated in S_CK, and related to purine metabolism, ubiquitin mediated proteolysis, and amino sugar and nucleotide sugar metabolism. Genes in cluster 3 mainly enriched in tyrosine metabolism, linoleic acid metabolism and alpha-linolenic acid metabolism, and showed lower expression in the S_CK group compared to other groups. For cluster 4, 19 genes were upregulated in the resistant cultivar (R_CK and R_T), but downregulated in the susceptible cultivar (S_CK and S_T). These genes were closely related to cutin, suberine and wax biosynthesis, sesquiterpenoid and triterpenoid biosynthesis, and steroid biosynthesis pathways. Genes in cluster 5 showed the opposite expression patterns to cluster 4. Speci cally, 20 DEGs were upregulated in Caoyuan No.2 (S_CK and S_T), but downregulated in the Caoyuan No.4 cultivar (R_CK and R_T). These genes encoded cysteine-rich receptorlike protein kinase (4 genes), rust resistance kinase and WRKY transcription factor. Cluster 6 contained 6 DEGs that were upregulated in the resistant cultivar (R_CK and R_T), encoding nodule lectin. However, most DEGs were unnamed, and had an uncertain function in the heat map, which needs further study to understand their functions (Fig. 5b, TableS4). The other 2087 constitutive defensive genes in Caoyuan No.4 were mainly enriched in energy metabolism (41 genes), carbohydrate metabolism (35 genes), replication and repair (125 genes), translation (35 genes), environmental adaption (22 genes), and lipid metabolism (16 genes) (Fig. 5c, Table S5).

Discussion
Alfalfa is an extremely energy e cient crop and is playing an increasingly important role in low input sustainable agriculture. However, more than 100 insect species damage alfalfa in southeast Asia, northeast Africa, and the U.S. [29]. Thus, the use of insect-resistant cultivars to control insects that damage both the quantity and quality of the alfalfa is an effective strategy. A wide range of population breeding methods including phenotypic recurrent selection have been especially successful in developing resistance to insects [30]. The high thrips-resistance alfalfa cultivar Caoyuan No.4 was also bred using recurrent selection in our previous study [21].
Plants have evolved effective defense mechanisms against insect infestation including mechanistic (trichomes, hairs) defenses and chemical defenses that involve genes and pathways related to diverse mechanisms [31]. For breeding insect resistance in crops, it is necessary to have information on genetic variation in the host reaction to insect infestation. In this study, we identi ed the gene functions and the molecular pathways related to thrips resistance through RNA-Seq analysis of the thrips-resistant alfalfa accession (Caoyuan No.4) and a thrips-susceptible alfalfa accession (Caoyuan No.2).
The results of PCA analysis con rmed that Caoyuan No.2 was more susceptible to thrips than Caoyuan No.4 (Fig. 2, 3). Consistent with this result, very few DEGs (only 13) existed in Caoyuan No.4 after thrips stress, but a total of 3326 contigs DEGs were detected in Caoyuan No.2 after thrips stress. Both results con rmed our previous results [21], indicating that Caoyuan No.4 was a thrips-resistant cultivar, while Caoyuan No.2 was a thrips-susceptible one. Thus, the choice of genotype and the reproducibility of the assay were suitable for subsequent DEGs analysis.
Every plant species has its own mechanism to protect itself from insect attack, including constitutive and induced defence mechanisms. In this study, 1285 constitutive differential genes were expressed between Caoyuan No.4 and Caoyuan No.2, while only 13 DEGs were induced by thrips attacks in Caoyuan No.4. In addition, the results of qRT-PCR and RNA-seq veri cation revealed that genes involved in cell wall biogenesis and metabolism expressed higher in thrips-resistant alfalfa accession (Caoyuan No.4) than those of thrips-susceptible alfalfa accession (Caoyuan No.2) (Fig. 6). These result con rmed our hypothesis that the molecular mechanisms of constitutive defence may be more important for thrips resistance in Caoyuan No.4.
In addition to the induced production of physical and chemical defenses, numerous changes in plant primary metabolism, such as carbohydrate and nitrogen metabolism, occur in response to insect attacks [32]. However, there are different theories as to how the plant primary metabolite responds to plant insect defense. It has been reported that plant genotypes that are able to maintain effective photosynthesis under insect attack often exhibit greater insect resistance [33,34], and increased photosynthesis or carbohydrate metabolism can serve as energy sources for the production of plant defenses in many plant-insect interactions [32]. Genes related to the energy and carbohydrate metabolism including photosynthesis (map00195), carbon xation (map00710) and nitrogen metabolism (map00910) were upregulated, indicating that the primary metabolism may play an important role in responding to thrips attack in Caoyuan No.4. Plant epicuticular lipid extracts and individual lipid components such as cutin and wax are important for plant insect resistance by affecting oviposition, movement, and feeding [35,36]. Many wax-related genes such as CER1, CER2 and CER3 have been identi ed in different plants [37,38]. In Caoyuan No.4, genes associated with lipid metabolism, including cutin, suberine and wax biosynthesis were also upregulated, suggesting that lipid metabolism may play an important role in the thrips-resistance of Caoyuan No.4. Analysis of leaf surface wax of Caoyuan No.4 and Caoyuan No.2 con rmed this hypothesis (unpublished data). Plant secondary metabolites such as alkaloids, glucosinolates or phenolic compounds serve as plant defenses to insects [39,40]. However, the content and distribution of individual secondary metabolites vary greatly among plant genotypes [41]. Here, the important roles of zeatin, sesquiterpenoid [42] and triterpenoid [43], carotenoid and diterpenoid [44], monoterpenoid [45] and phenylpropanoid [46] were also highlighted by the existence of a large number of upregulated genes involved in the synthesis of these chemicals in Caoyuan No.4, which indicated that thrips defense in Caoyuan No.4 was related to the existence of certain secondary metabolites. It is well known that the interaction of plant-pathogen and plant-insects share some strategies. Although the genes related to plant-pathogen interaction (CDPK, CaMCML, Pit 6 , RPM 1 , CERK1, SGT, HSP 90 , EDS 1 , WRKY, NHO 1 and PR) [47] were signi cantly upregulated in Caoyuan No.2 after thrips stress, the genes (CDPK, CNGC S , FLS 2 , RPM 1 and EDS 1 ) [48] involved in plant-pathogen interaction were signi cantly downregulated in Caoyuan No.4 compared to Caoyuan No.2 (Sup.S1). Comparing the genes related to plant-pathogen interaction expressed in the two groups, we found that the FLS 2 gene only expressed in Caoyuan No.4, which is related to cell wall reinforcement by ROS reaction. Moreover, more genes related to plantpathogen interaction expressed in Caoyuan No.2 after thrips infestation. Consistent with our previous studies and other reports, the results indicated that Caoyuan No.4 was resistant to thrips because of reinforcement of the cell wall [20], while Caoyuan No.2 could develop thrips defense by inducing plantpathogen interaction genes such as NHO1and PR1 [49]. In addition, genes related to iso avonoid and avonoid biosynthesis, proteasome, amino sugar and nucleotide sugar metabolism, carbohydrate metabolism and energy metabolism were all upregulated in Caoyuan No.2 after thrips infestation, indicating that these genes might be induced defences in responding to thrips attack. These ndings indicated that the defence mechanism in Caoyuan No.4 was mainly constitutive resistance, including primary metabolism, lipid metabolism and certain secondary metabolites, but not the plant-pathogen interaction associated genes.

Conclusion
Both constitutive defensive genes and potential induced defensive genes were detected in the defense of Caoyuan No.4 against thrips attack. This nding is signi cant for plant resistance breeders, because the two distinct kinds of defensive genes -constitutive defensive genes and induced defensive genes -can be simultaneously activated and thus potentially enhance plant protection against insects attacks.