Full-length transcriptome and targeted metabolome analyses provide insights into defense mechanisms of Malus sieversii against Agrilus mali

Plant materials

On June 25, 2017, 50 shoots from 10 selected trees (five resistant and five susceptible) were collected from the wild M. sieversii population located in Gailiangchang (N43°22, E83°34, the altitude 1,280–1,450 m), Xinyuan County, Xinjiang Autonomous Region. Field experiments  was approved by the research council of the Institute of Horticultural Crops, Xinjiang Academy of Agricultural Sciences (approval number: XAASHC2016.02.08). The selected trees (resistant and susceptible) have been proved to be stable in their insect-interaction traits for years. All the trees were aged 45–60 years old.

Infection treatment

The apple buprestid beetles usually parasitize the phloem of branches. According to their life cycle, adults lay eggs in the phloem underneath the bark. The larvae begin feeding the follow spring. The damage caused by larval feeding peaks in July and lasts until September. Ten mature larvae (one larvae for per shoot) were inserted into each tree through a simulated pore canal that created by T shape cutting of phloem. After a period of ten days, phloem were sampled from the infested and uninfested parts of the resistant and susceptible trees after the  areas from the parts >3 cm away from the original mechanical damage had healed. (The sampled spots with obvious infection were chosen). Samples were frozen immediately in liquid nitrogen and stored at −80 °C before RNA extraction. Samples were named as follows: Insect-feeding resistant plant (IFR): Insect-feeding susceptible plants (IFS); None-feeding resistant plants (NFR); None-feeding susceptible plants (NFS).

RNA extraction, cDNA library construction, and Illumina sequencing

Total RNA was extracted from each sample using a RNA extraction kit (TissueLyzer; Qiagen, Valencia, CA, USA), according to manufacturer’s instructions. RNA purity was checked using a Nano Photometer spectrophotometer (IMPLEN, CA, USA), degradation and contamination were monitored on a 1% agarosegel, and integrity was assessed by an Agilent Bioanalyzer Nano 6000. Three µg of RNA was taken from each qualified sample and mixed. mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. The sequencing library was prepared with a NEBNext®Ultra™ RNA Library Preparation Kit (NEB, USA). First strand cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-). DNA polymerase I and RNase H were added sequentially into the obtained first strand cDNA to produce the second cDNA strand. The fragments in the library were purified using an AMPureXP system (Beckman Coulter, Beverly, USA). The obtained double-stranded cDNA was subjected to end repair, addition of an A-tail, and adaptor ligation, then AMPure XP beads were used for size selection and, finally, PCR was adopted for enrichment to construct a cDNA library. Library quality was assessed with an Agilent Bioanalyzer 2100 system. After qualification and quantification of the library, short-read sequencing of each treatment × phenotype combination was performed using the Illumina HiSeqX-ten platform for RNA-Seq analysis.

Construction of single molecule real-time library and sequencing

PacBio RS II and Illumina sequencing were done on the same sample from the same site. We mixed equal amounts of RNA and sequenced them. A SMARTer™ PCR cDNA Synthesis Kit was used to get full-length cDNA of mRNA, which was then screened with BluePippin™ Size. Three cDNA libraries (1–2 k, 2–3 k, 3–6 k) were constructed. PCR was used again to amplify the selected full-length cDNA. Then, they were repaired at the end, ligated with the SMRT dumbbell-type adapters, and cleaved by exonuclease. A secondary screening was performed using Blue Pippin to obtain the sequencing libraries. After the libraries were constructed, they would quantify accurately using Qubit2.0, and the fragment size of libraries was detected using Agilent 2100. Only those complying with a standard could be sequenced. Long-read sequencing was performed using the PacBio RS II platform.

Preprocessing of SMRT reads

Reads of Insert (ROI) sequences were extracted using the following criteria: > = 0 full-pass reading of the polymerase and sequence accuracy > = 75. Readings below 50 bp in length were discarded. ROI sequences were divided into two categories: full length sequences (that included a 5′ primer, a 3′ primer, and a polyA tail) and non-full-length sequences. We used the ICE algorithm of SMRT Analysis (v2.3.0) for iterative clustering, which collected similar sequences (i.e., multiple copies of the same transcript) into the same cluster. The parameters used are the default parameters. We used quiver to correct consensus sequences of each cluster, which yielded high quality transcripts (HQ) with an accuracy >99% (Korlach et al., 2010). In this project, the low-quality consensus sequences obtained from each sample were corrected using the corresponding Illumina RNA seq data through proof read to improve sequence accuracy (Hackl et al., 2014). PacBio RS II and Illumina sequencing data are deposited at the National Center for Biotechnology Information (NCBI) database (NCBI Sequence Read Archive SRP199016).

Function annotation and classification

The annotation of transcriptome sequences were performed using the following seven databases:  Nr (NCBI non-redundant protein sequences), Pfam (Protein family), KOG/COG (EuKaryotic Ortholog Groups and Clusters of Orthologous  Groups of proteins), Swiss-Prot (a manually annotated and reviewed protein sequence database), KEGG (Kyoto Encyclopedia of Genes and Genomes)  and GO (Gene Ontology). The read numbers were transformed to FPKM (fragments per kilobase of transcript sequence  per millions base pairs sequenced) value for gene expression quantification.

Real-time quantitative PCR verification

We applied qRT-PCR technology to validate the results of the differential expression analysis. Insect-resistant related genes with stable and drastic expression differences were selected for expression verification. Total RNA from each sample was reverse-transcribed using a 5× All-In-One MasterMix (Cat# G492, Abm, Canada) kit to synthesize the first strand of cDNA. The expression data in the Illumina HiSeq sequencing results were extracted for comparative analysis. Specific primers were designed using Primer-BLAST online software, and the melting curve was drawn to confirm PCR specificity. The quantitative PCR reaction conditions were as follows: 94 °C for 60 s, followed by 45 cycles of 94 °C for 5 s, 60 °C for 15 s, and 72°C for 10 s. SYBR was adapted to premix EvaGreen Express 2× qPCR MasterMix (Cat# MasterMix-ES, Abm, Canada). The test was performed with apple MdTUB1 as the reference gene using LightCycler 96 (Roche, USA). Three biological replicates were set for each group, and we used three technical replicates for each biological replicate. Please see Table S1 for the sequence of qRT-PCR primers.

Detection of secondary metabolites

The plant material was stored at −80 °C and ground with a grinder(MM 400, Retsch) at 30 Hz for 1 min. 100 mg of the powder was extracted in 1.0 mL of methanol, during which the sample was rotated three times to mix the solvent. After the extraction, 12,000g was centrifuged for 15 min to obtain a supernatant. The obtained supernatant was blown dry with nitrogen gas in a 35°C heating mode, then, resuspended in a mixture of 30% methanol with vortexing, centrifuged at 12,000 g for 15 min, and the supernatant was collected and stored for analysis on an LC/MS. The content of secondary metabolites was detected by liquid chromatography-mass spectrometry with reference to the method of Francescato et al. (2013), Pan, Ruth & Wang (2010) and Wojakowska et al. (2013). In this experiment, all chemical reagents were of analytically pure or chromatographically pure preparations. Chemical standards are from BioBioPha (http://www.biobiopha.com/) and Sigma Aldrich (http://www.sigmaaldrich.com/united-states.html).

The analysis conditions were as follows: The sample extracts were analyzed using an  LC-ESI-MS/MS system (HPLC, Shim-pack UFLC Shimadzu CBM30A,  http://www.shimadzu.com.cn/; MS,  Applied Biosystems 4500 Q TRAP,  http://www.appliedbiosystems.com/). The analytical conditions were as follows:  HPLC: column, Waters ACQUITY UPLC HSS T3 C18 (1.8 µm, 2.1 mm*100 mm); solvent system, water (0.04% acetic acid):  acetonitrile (0.04% acetic acid); gradient program, 95:5 V/V at 0 min, 5:95 V/V at 11.0 min, 5:95 V/V at 12.0 min, 95:5 V/V at 12.1 min,  95:5 V/V at 15.0 min; flow rate, 0.40 ml/min; temperature, 40 °C; injection  volume: 5 µl.  The effluent was alternatively connected to an ESI-triple  quadrupole-linear ion trap (Q TRAP)-MS.

For quantification of metabolites, standard solutions of different concentrations were prepared. Then mass spectrometry was used to analyze 0.1 ug mL⁻¹, 0.2 ug mL⁻¹, 0.5 ug mL⁻¹, 1.0 ug mL⁻¹, 2.0 ug mL⁻¹, 5.0 ug mL⁻¹, 10.0 ug mL⁻¹standard solutions and samples.

Correlation analysis of transcriptome and metabolome data

By combining the metabolomics and transcriptomics data through the KEGG metabolic pathway, differentially expressed/accumulated genes and metabolites involved in the same biological process (KEGG Pathway) can be identified quickly.

Statistical analysis

All the experimental data obtained were analyzed using SPSS 17 software (SPSS, Inc.  Chicago, 186 IL, USA) and indicated by means  ± (n = 3) standard deviation (SD). Data  were analyzed using One-way ANOVA and Duncan test at a significance level of  p < 0.05.

Morphological and physiological differences between resistant and susceptible plants

Susceptibility to the Agrilus mali invasion, soluble sugar, polyphenol as well as tannin contents of M. sieversii is shown (Fig. 1) to illustrate the difference between resistant and susceptible plants. Greater contents for polyphenol and tannin were observed in resistant plant shoots compared to susceptible plants (Figs. 1A–1B). However infection rate (number of insects) and soluble sugar contents in the shoots of susceptible plants were significantly higher than resistant plants (Fig. 1D). As the important secondary metabolites in plants, polyphenols and tannin are closely related to the disease resistance of plants, also play an important role in plant-insect interaction. Therefore higher polyphenol and tannin contents in resistant plant shoots may be related to resistance of these plants to Agrilus mali invasion. In the later section, we will focus on pathways involved biosynthesis of secondary metabolites such phenolic acids and plant-pathogen interaction.

PacBio RS II and Illumina HiSeq sequencing

Transcriptome sequencing was performed using Illumina HiSeq X ten and PacBioRS II. First, all RNA samples were tested with an Illumina HiSeq X ten, which resulted in a total of 130.06 Gb Clean Data and Q30 up to 93.38% (Table S2). Then, an equal amount of the 12 RNA samples were mixed together for PacBio library preparation and sequencing. The full-length cDNA of the poly(A) + RNA sample was standardized and SMRT-sequenced using the PacBio RS II platform. The number of sequencing cells was two for a 1–2 K library, and one each for the 2–3 K and 3–6 K libraries. A total of 122,103 full-length non-chimeric (FLNC) reads (Table 1) and 275, 586 ROI (reads of insert) (Table 2) were obtained. SMRT Analysis software was employed to cluster the full-length sequences, and a total of 64,692 consensus transcript sequences were generated. Non-full-length sequence correction produced 47,181 high-quality, full-length transcript sequences and 17,511 low-quality sequences, which were corrected using Illumina Hiseq data (Table S3).

The sequences were clustered iteratively using SMRT Analysis (v2.3.0) software and the ICE (Iterative Clustering for Error Correction) algorithm. In this experiment, we received a total of 64,692 F01 consensus sequences with HQ (High-quality) and LQ (low-quality) transcripts that were obtained from different libraries (Table S3). At least 27,230 HQ transcripts were obtained from the 1–2 kb library, and 9,895 HQ were obtained from the 2–3 Kb library. In addition, the LQ consensus sequences were corrected by Illumina RNA seq platform using proof reads software to improve their accuracy (Hackl et al., 2014).

Verification of RNA-Seq gene expression

To confirm the reliability of the RNA-Seq data, synthetic pathways, transcription factors, and important genes related to plant defense were screened from the annotation library of RNA-Seq DEGs data, and 18 of them were selected for real-time quantitative PCR (qRT-PCR). By doing so, effects of infection on expression of defense-related genes, and the accuracy of RNA-Seq data were verified (Fig. 2). Quantitative results showed that all 18 genes responded to infection intensively, which included 13 up-regulated and five down-regulated genes. Most changes in gene expression were consistent with the expression patterns in RNA-Seq results, which confirmed the reliability of Illumina HiSeq sequencing in this experiment. Furthermore, expression level of a gene annotated as protein TIFY9 increased significantly after insect induction, which exhibited more than a 30-fold up-regulation in resistant strains and a more than 10-fold up-regulation in the control. This indicated the active participation of TIFY genes in the regulation of an insect-resistant response.

Functional annotation and selection of candidate DEGs related to plant pest resistance

In order to identify genes respond to Agrilus mali invasion in resistant plants, firstly we select common DEGs in NFS vs NFR (A), IFS vs IFR (B), NFR vs IFR(C) and NFS vs IFS (D) comparison groups, respectively. Then we used these common DEGs to build the Venn-diagram (Fig. S2). It is evident from the Venn graph that 266 DEGs respond in resistant genotype compared to the susceptible one and the control (non-feeding) samples(Fig. 3A). These DEGs mainly enriched in the pathways of phenylalanine metabolism, degradation of ammonic acids, Flavonoid biosynthesis, phenylpropanoid biosynthesis, terpenoid backbone biosynthesis, degradation of aromatic compounds, fatty acid metabolism, phosphatidylinositol signaling, plant hormone signal transduction pathway, plant-pathogen interaction etc. (Fig. 3B). Among these DEGs 240 of 266 were up-regulated and 26 DEGs were down-regulated in resistant plants after infection by Agrilus mali indicating that these genes may participated in plant-insect interaction (Fig. S3). In addition majority of these genes enrichen in secondary metabolic pathways verifying the important role of secondary metabolites in plant defense system. Among these 266 DEGs, we mined 19 DEGs that enriched in phenylalanine metabolic pathway, flavonoid biosynthetic pathway and plant-pathogen interaction. These selected 19 genes will be further analyzed in following section.

Correlation analysis of transcriptome and metabolome

Ten days after insect feeding, we analyzed a group of bioactive secondary metabolites such as phenolic acids, alkaloids, flavones, salicylic acid, chlorogenic acid and benzoic acid that involved in plant defense system (Fig. S4). Among these metabolites, only 3 phenolic acids differentially accumulated in different samples (Figs. 4A–4C). Trans-cinnamic acid and caffeic acid content increased with 1.3-fold, 1.1-fold and 7.3-fold, 2.8-fold in susceptible and resistant plants, respectively. The fold changes of these two metabolites were much higher in resistant plant than in susceptible plant, indicating that these metabolites may play critical role in plant-insect interaction. Through correlation analysis, we found that genes encoding enzymes (PAL, C4H) responsible for biosynthesis of these four phenolic acids were differentially express in susceptible and resistant plants (Fig. 4D). In accordance with metabolites content, expressions of correlated genes were much higher in resistant plants compared to susceptible plants.