Metabolomic and transcriptomic analyses of Prunus davidiana reveal the key resistance compound to green peach aphid, betulin

Prunus davidiana, a close wild relative of the cultivated peach, has been identied as having a strong resistance to Myzus persicae Sülzer (green peach aphid, GPA), one of the major pests of peach. However, the resistance mechanism of P. davidiana remains unclear. In this study, combined analysis of metabolome and transcriptome was conducted using aphid-resistant (R-32) and aphid-susceptible (S-27) lines, from a segregating population, to investigate the defense mechanism of P. davidiana. These results showed that R-32 continuously accumulated higher levels of betulin than S-27 during aphid infestation. Besides, betulin displayed a strong toxic effect on GPA. Transcriptome analysis revealed that the expression of genes involved in the betulin biosynthesis pathway responded highly to GPA infestation, especially CYP716A1. Our results demonstrate that betulin play an important role in the mechanism of resistance of P. davidiana to GPA.


Introduction
Crop wild relatives are abundant germplasm resources in nature and possess bene cial genes that are missing from their cultivated crops, especially resistance genes. The study of wild relative resistant germplasms with excellent agronomic traits can improve the understanding of the resistance of crops to environmental stress and promote the development of sustainable agriculture by reducing the use of chemicals 1 . Peach (Prunus persica (L.) Batsch) constantly encounters various unfavorable environments during its growth, including abiotic and biotic stresses. The green peach aphid (GPA, Myzus persicae Sülzer) is one of the major pests of peach, which not only affects its growth and production, but also mediates the spreading of a series of plant viruses. However, the gene that provides resistance to GPA is absent in the majority of currently cultivated peach varieties. Several aphid-resistant peach genotypes have been identi ed in previous studies. In the last century, two genotypes have been identi ed, namely "Weeping Flower Peach" (P. persica) and "Rubira" (P. davidiana) [2][3][4] . Recently, Niu et al. have shown that "Fen Shouxing" (P. persica var. densa Makino) was also an aphid-resistant genotype and speculated that the WRKY DNA-binding protein gene (Prupe.1G564300) may be the gene that provides resistance 5,6 . However, the defense mechanism of P. davidiana, the variety with the strongest resistance among the 15 aphid-resistant germplasms in China 7 , is not clear. P. davidiana, as a close wild relative of the cultivated peach, has been reported to have an antibiosis effect on GPA, which means that the fecundity of GPA decreased and the mean length of the mature embryos within the gonads of the females on the day of adult molt was negatively correlated with the total number of embryos 8 . Therefore, there are certain metabolites of P. davidiana that confer resistance to GPA. However, it is still not clear which metabolites of P. davidiana play a role in GPA resistance. The recently developed metabolomic technology has allowed the detection of metabolites at a highthroughput level. Several studies have identi ed metabolite changes during the infestation of Ostrinia furnacalis, Mythimna separata, Macrosiphoniella tanacetaria and Uroleucon tanaceti on maize and Tanacetum vulgare, respectively, using metabolomic technology [9][10][11] . However, the metabolomic changes during peach infestation by aphids remain unclear. Sauge 12 showed that the major aphid-resistance quantitative trait loci (QTL, MP.SD-3.1), on G3 close to the AG106 and AG50B RFLP markers, was identi ed in P. davidiana. However, the key gene for GPA resistance in P. davidiana remains unclear.
Zhouxingshantao 1, a representative germplasm of P. davidiana, shows strong resistance to aphids. In the present study, D29-32 (R-32) and D29-27 (S-27) from the segregating population of P. davidiana var. Zhouxingshantao 1 were resistant and susceptible, respectively, and were used to investigate the defense mechanism of P. davidiana. To unravel the underlying molecular defense mechanisms, we performed an integrated metabolomic and transcriptomic analysis of both R-32 and S-27 during aphid infestation to identify candidate key metabolites and genes that contribute to defense responses. This work provides a solid foundation for further research on the major genes responsible for resistance to GPA in P. davidiana.

Results
The rst two days post-infestation were the key period for the GPA response in P. davidiana Apterous GPAs were settled on the partially expanded leaves of ve annual new shoots of resistant (R-32) and susceptible (S-27) individuals two weeks after spraying the chemical insecticide. The settlement and number of aphids per shoot were recorded every day for one week after aphid infestation (Fig. 1a, b). The number of aphids on R-32 slowly decreased and almost disappeared from 4 days post-infestation (DPI) and after. However, the aphid population on S-27 increased rapidly and gradually stabilized from 2 DPI and after. At the same time point, the number of aphids on S-27 was signi cantly higher than that on R-32 (p < 0.01). Therefore, we hypothesize that the resistance of R-32 to GPA is due to changes in metabolites in the rst two days post infestation.
The differentially accumulated metabolites (DAMs) (Log 2 FC (fold change) ≥ 1 and VIP ≥ 1) during GPA infestation (from 0 to 48 hpi) in resistant and susceptible accessions were identi ed by the comparison between GPA treatments and control (Fig. 1c). Finally, 192 and 217 DAMs were identi ed during GPA infestation in R-32 and S-27, respectively. Interestingly, the down-regulated DAMs were much more abundant than the up-regulated DAMs in both R-32 and S-27 in each treatment, suggesting that metabolites may participate in the response to GPA infestation via negative regulation. Notably, the changes in metabolite expression peaked at 6 hpi in R-32 and at 12 hpi in S-27, suggesting that the responses to aphid infestation in resistant accessions were faster than in susceptible accessions.
To study the changing trends of metabolite relative contents during GPA infestation, K-means cluster analysis was performed using Pearson correlation distances based on the expression pattern during GPA infestation. Finally, the 341 metabolites, with differences in at least two time points, could be classi ed into six major groups with 12 subclasses (Supplementary Fig. 1). Metabolites in group 1 (including 38 metabolites) was dominated by the avonoids (~ 47.37%), harboring stable expression throughout the aphid infestation with the expression level in R-32 signi cantly higher than that in S-27. Group 2 could be divided into two subclasses, subclass 2 (34 metabolites) and subclass 3 (18 metabolites). Intriguingly, the expression of these metabolites peaked at 6 hpi in R-32, while at 12 hpi in S-27. The expression of metabolites in group 3, including subclass 4 (52 metabolites) and 5 (15 metabolites), showed different pattern in R-32 and S-27 lines, with a "delince-increase-decline" pattern in R-32 line and a "decline" pattern in S-27 line. In group 4, including subclass 6 (30 metabolites) and subclass 7 (19 metabolites), these metabolites were up-regulated at both 6 hpi and 24 hpi in R-32 and up-regulated at 12 hpi in S-27. The metabolites in group 5 (76 metabolites) expressed stably in R-32 line but changed dramatically in S-27 line. In group 6, including three subclasses 10 (35 metabolites), 11 (9 metabolites), and 12 (15 metabolites), the expressions of these metabolites in R-32 were down-regulated at early stage of GAP infestation but up-regulated at late stage. For S-27, the lowest expression level was at 12 hpi for subclass 10 and 11, while the continuous increase was observed in subclass 12. Differences in metabolite expression during GPA infestation may be part of the defense mechanisms of P. davidianan var. Zhouxingshantao 1. From all six groups, it is interesting to observe that metabolites in Group 1 had higher expression levels in R-32 than in S-27 throughout the infestation by GAPs. These may be inherent differences between the two accessions, which include a category of metabolites that contribute to the constitutive defense mechanism. Candidate causative metabolites involved in GPA resistance To identify the possible causative metabolites for GPA resistance in peach, we compared the DAMs between R-32 and S-27 at all ve time points (Fig. 2). In total, we identi ed 59, 81, 50, 68, and 61 upregulated DAMs in the comparison of R-32 and S-27 at 0, 6, 12, 24, and 48 hpi, respectively. Among these, 21 were shared between all the comparisons. Among the 21-shared DAMs, 19 metabolites belonged to subclass 1, and two metabolites belonged to subclass 4. Similarly, 152 downregulated DAMs were identi ed in the comparison of R-32 and S-27 during GPA infestation. Of these, 11 were shared at the ve time points, including eight metabolites from subclass 9 and three metabolites from subclass 8 (Fig. 2a).
To further lter the representative metabolites in the GPA-resistant mechanism, we compared the top 20 up-and down-regulated DAMs between R-32 and S-27 at all ve time points (Fig. 2b-2f). We found that betulin had the highest Log 2 FC (> 15) throughout the experiment, followed by quercetin 7-Omalonylhexosyl-hexoside (Log 2 FC > 13). Such high Log 2 FC may be due to the very high levels of these two metabolites in R-32 compared to S-27. Although studies on the anti-insect effects of betulin are rare, betulin derivatives have been shown to have a strong effect against insects. Betulinic acid, an insect growth regulator, has been shown to effectively interfere with the growth of Spodoptera litura 16 and Callosobruchus chinensis 17 . Therefore, it is reasonable to hypothesize that betulin, a lupane-type triterpene, may have an effect on GPAs. Quercetin, one of the most abundant avonoids in plants, has been documented to be important in plant-insect interactions [18][19][20] . All these studies support our hypothesis that the quercetin derivative, quercetin 7-O-malonylhexosyl-hexoside, plays a role in the resistance to GPAs in P. davidiana. High toxicity of betulin to the GPA As glycosylation and malonylation of avonoids increase their water solubility and stability 21,22 we used quercetin instead of quercetin 7-O-malonylhexosyl-hexoside to verify its insecticidal toxicity. The toxicity of betulin and quercetin to GPA was determined using the leaf dipping method. The median lethal concentration (LC 50 ) values of betulin and quercetin after 12 and 24 h of exposure were further calculated (Fig. 3a). The results showed that LC 50 values of betulin and quercetin after 24 h of exposure were 396.08 and 727.56 mg⋅L − 1 , respectively. After 48 h exposure, the LC 50 values of betulin and quercetin were 206.64 and 306.34 mg⋅L − 1 , respectively. This implies that betulin has a stronger toxic effect on GPA than quercetin. To further evaluate the inhibitory effect of betulin against aphids, we recorded the number of aphids on S-27 shoots after spraying with acetone (S-27 + acetone) or betulin (396.08 mg⋅L − 1 , S-27 + betulin) (Fig. 3b). This experiment showed a similar slow increase in the number of aphids on S-27 + acetone to that on S-27. Compared to S-27 + acetone, the number of aphids on S-27 + betulin remained at a lower level. After aphid infestation for 7 days, dead aphids were observed on S-27 + betulin (Fig. 3d). This indicated that betulin had an inhibitory effect on aphids. Moreover, the absolute content of betulin in R-32 (66.59 µg/g DW) was signi cantly higher than that in S-27 (55.24 µg/g DW) (Fig. 3c). Therefore, we suggest that betulin may play a crucial role in the resistance of P. davidiana to GPA.
Transcriptomic changes during the GPA infestation in peach To elucidate the GPA defense mechanisms at the transcription level, RNA-seq was performed for R-32 and S-27 using the same samples used for the metabolite analyses. We generated 916 million clean reads using the Illumina platform. After aligning against the peach reference genome "Lovell" (release version 2.0) 23 , the average mapping rate was calculated to be 95.40% (Supplementary Table 2). Finally, 17, 813 genes expressed (FPKM 1) in at least one sample were identi ed, of which 13, 923 were expressed in all samples (Supplementary Table 3). The RNA-seq results were further veri ed by qRT-PCR ( Supplementary   Fig. 2). A total of 5, 563 differentially expressed genes (DEGs, |Log 2 (FC)| > 1) were identi ed, of which 1, 945 and 1, 824 were differentially expressed during GPA infestation in R-32 and S-27, respectively. In addition, 4,378 genes showed differential expression in the comparison of R-32 and S-27 during the whole experiment, suggesting the presence of differences at the transcriptional level between resistant and susceptible accessions.
DEGs during GPA infestation (from 0 to 48 hpi) in resistant and susceptible accessions were identi ed by comparing GPA treatments and the control (Fig. 4a). Notably, the DEGs were abundant at 6 hpi and 12 hpi for both R-32 and S-27. To further identify which known biological pathways were enriched with DEGs during GPA infestation, we performed KEGG pathway enrichment analysis with DEGs between the resistant and susceptible accessions in the rst two days post aphid infestation (Fig. 4b-4f). Without aphid infestation, DEGs involved in 6 pathways were enriched. Among these pathways, the fold enrichment of avonoid biosynthesis pathway was the highest. There were 21 pathways enriched with DEGs at 6 hpi. Linoleic acid metabolism pathway was highly enriched, second by brassinosteroid biosynthesis and monoterpenoid biosynthesis pathways. At 12 hpi (15 pathways enriched), avonoid biosynthesis pathway was signi cantly enriched. Seven pathways were enriched with DEGs at 24 hpi. Phenylpropanoid biosynthesis had the highest fold enrichment, followed by avonoid biosynthesis pathway. At 48 hpi, there were 14 pathways were enriched. DEGs involved in stilbenoid, diarylheptanoid and gingerol biosynthesis pathway was highly enriched, second by avonoid biosynthesis pathway.
Except to biosynthesis of secondary metabolites pathway, avonoid biosynthesis and phenylpropanoid biosynthesis pathways were continuously enriched by DEGs during aphid infestation. As phenylpropanoids are important precursors for the synthesis of avonoids, this result was consistent with the accumulation of quercetin derivatives. However, the DEGs were not enriched in the terpene biosynthesis pathway. This may be due to the complexity of the terpenoid biosynthesis pathway and the unannotated biosynthesis pathway of betulin in KEGG. Therefore, the related biosynthesis pathway of betulin was not found to be enriched.

Expression pattern of genes in betulin pathway during GPA infestation
To elucidate the GPA defense mechanisms at the gene level, we investigated the expression of genes involved in the betulin biosynthesis pathway. However, the betulin biosynthesis pathway in peaches remains unclear. Based on the production process of betulin in transgenic yeast 24 , we selected possible genes involved in betulin biosynthesis to investigate the expression pattern of the genes at 0, 6, 12, 24, and 48 hpi in both R-32 and S-27. A total of 19 genes were selected, which were expressed (FPKM > 2) in the two genotypes at least at one time point (Fig. 5a). The relative expression of the 19 genes was normalized using logarithmic standardization. This showed that the expression of genes involved in the betulin biosynthesis pathway responded highly to GPA infestation. HMGCS (Prupe.5G088900), HMGCR2 (Prupe.8G182300), FDPS1 (Prupe.4G002700), FDPS2 (Prupe.6G028100), FDFT1 (Prupe.8G08700), SQLE3 (Prupe.8G156800), LUS1 (Prupe.3G025900), CYP716A1 (Prupe.1G002500), CYP716A3 (Prupe.4G103000) and CYP716A5 (Prupe.4G103200) were all up-regulated at 12 hpi for both R-32 and S-27. All of these genes, apart from HMGCR2 and FDPS1, were expressed at higher levels in R-32 than in S-27 at 12 hpi. In particular, the expression of CYP716As in R-32 was consistently higher than that in S-27 during aphid infestation. In addition, the FPKMs of these CYP716As were more than twice in R-32 compared to S-27 at 6 hpi and 12 hpi, which was veri ed by qRT-PCR ( Supplementary Fig. 2). This indicated that CYP716As might have an effect on the resistance of R32 to GPA.

Discussion
There are detailed studies on aphid resistance in maize 28 , chrysanthemum 29,30 , pepper 31 , and others, however, studies on the defense mechanism of P. davidiana to aphids are relatively lagging. Although the genetics of P. davidiana resistance to aphids have been studied since the 1990s 8 , research on the underlying molecular resistance mechanism lags far behind. This study describes the rst effort in the comparative analysis of metabolomics and transcriptomics to identify the genes and metabolites involved in P. davidiana responses to the infestation of Myzus persicae Sülzer. The results shed light on the possible mechanisms of the dynamic responses of P. davidiana to Myzus persicae Sülzer.
Betulin may have an effect on the resistance of P. davidiana to GPA In this study, metabolomic analysis revealed that the relative concentrations of betulin and acylglucosidate quercetin in R-32 were continuously higher than that in S-27 during aphid infestation. However, the terpene biosynthesis pathway was not enriched during aphid infestation. This may be due to the complexity of the terpenoid biosynthesis pathway and the unannotated biosynthesis pathway of betulin in KEGG. Furthermore, bioassays and eld validation in the greenhouse showed that betulin had a stronger inhibitory effect on aphids than quercetin. It is reasonable to hypothesize that betulin may play a crucial role in the resistance of P. davidiana to GPA, while acylglucosidate quercetin may play a supporting role.
Betulin, a lupane-type triterpene, has received wide attention for its broad spectrum of biological activities. It has been demonstrated that betulin and betulinic acid possess antimalarial, antiin ammatory, antibacterial, antitumor, and anti-HIV properties [32][33][34][35][36] . However, studies on the anti-insect effects of betulin are rare. At the end of the last century, several studies have reported that betulin derivatives have a strong antifeedant effect against bollworm larvae (Heliothis zea) and Colorado potato beetle (Leptinotarsa decemlineata), and betulinic acid derivatives against tobacco caterpillars (Spodoptera litura F) 37 . Recently, betulinic acid, an insect growth regulator, has shown effective interference with the growth of Spodoptera litu 16 and Callosobruchus chinensis 17 . This suggests that betulin may play a crucial role in the resistance of P. davidiana to GPA.

CYP716As may play a role in regulating betulin biosynthesis
As the identi cation of plant genes involved in betulin synthesis, its biosynthesis pathway in plants is becoming increasingly clear. Recent studies have revealed that betulin can be produced in yeast by catalyzing lupeol with CYP716A, a subfamily of cytochrome P450 monooxygenases 24,38,39 . In this study, we identi ed homologs of CYP716A175 24 in peach and selected three expressed CYP716As (FPKM > 1, at least one time point), including CYP716A1, CYP716A3, and CYP716A5. The expression patterns of genes involved in the betulin biosynthesis pathway indicated that this pathway highly responded to GPA infestation, and CYP716As may be the key genes that mediate betulin biosynthesis.
As the largest enzyme family in plant metabolism, cytochrome P450 plays a crucial role in the detoxi cation of xenobiotics and the biosynthesis of secondary metabolites, antioxidants, and phytohormones in higher plants 40 . Among the cytochrome P450 superfamily, the CYP716A subfamily is mainly involved in triterpene biosynthesis 41,42 . Recently, Medicago truncatula CYP716A12 was identi ed as a multifunctional enzyme with the ability to catalyze consecutive oxidation at the C-28 position of αamyrin, β-amyrin, and lupeol to produce ursolic, oleanolic, and betulinic acids, respectively 24,26 . Its homologs, Vitis vinifera CYP716A15 and CYP716A17, catalyze the oxidation of lupeol to betulin and betulinic acid 24 . Lotus japonicus CYP716A51 also exhibits triterpenoid C-28 oxidation activity and participates in betulinic acid biosynthesis 27 . Phylogenetic analysis revealed that P. persica CYP716A1 is clustered into a clade with Medicago truncatula CYP716A12 and Lotus japonicus CYP716A51. This implies that CYP716A1 may have a function similar to Medicago truncatula CYP716A12 and Lotus japonicus CYP716A51. Therefore, CYP716A1 may play a role in regulating betulin biosynthesis.
To investigate the molecular mechanisms of resistance to GPA in P. davidiana, integrated analyses of transcriptomic and metabolomic pro les of the resistant line R-32 and the susceptible line S-27 during the infestation of GPA were performed. The results showed a different metabolomic pattern between the resistant and sensitive peach accessions and revealed a strong toxic effect of betulin on GPA. Therefore, we proposed a model of P. davidiana resistance to green peach aphid (Fig. 6). Analyses of the expression of genes involved in betulin biosynthesis identi ed CYP716As as strong candidate genes regulating the expression of betulin. This study provides new insights into the prevention of aphid infestation and candidate key genes involved in resistance to GPA in P. davidiana.

Materials And Methods
Plant materials, aphid infestation, and sampling Two ve-year-old accessions from an F 2 progeny derived from a cross between P. persica var.
For aphid infestation, treatment was performed at two weeks (pesticide duration) after spraying the insecticide Closer (Dow AgroSciences, USA). A total of 15 adult GPAs (Myzus persicae Sülzer), collected from a naturally infested peach orchard (Zhengzhou Fruit Research Institute, Zhengzhou, China), were placed on the partially expanded trifoliate leaves of each shoot for each individual and subsequently covered with an argenteus y net. Shoot tips with young leaves were sampled at time 0 (no aphid infestation), 6, 12, 24, and 48 h post infestation (hpi). Before sampling, the remaining aphids were carefully removed from the leaves. The collected shoot materials were immediately frozen in liquid nitrogen and stored at -80°C. Two biological replicates were used for each sampling point. RNA isolation, mRNA sequencing, and RNA-Seq Data Analysis The same samples used for the metabolome analyses were also used to perform RNA-sEq. Total RNA was extracted using the Quick RNA Isolation Kit (Huayueyang Biotechnology Co., Ltd., Beijing, China) according to the manufacturer's instructions. RNA purity and RNA integrity were examined using the kaiaoK5500®Spectrophotometer (Kaiao, Beijing, China), and the concentration was con rmed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). A total of 24 sequencing libraries (two biological replicates at six time points for both the aphid-resistant andsusceptible lines) were created from 2 µg RNA per sample using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (#E7530L, NEB, USA) following the manufacturer's recommendations. RNA concentration of libraries was measured using Qubit® RNA Assay Kit in Qubit® 3.0, and each library's quality was then assessed using the Agilent Bioanalyzer 2100 system (Agilent Technologies).
Reference peach genome v2.0 and annotation les were downloaded from the GDR website ( https://www.rosaceae.org/species/prunus_persica/genome_v2.0.a1) 43 . The genome index was built using Bowtie2 v2.2.3 44 . Clean data (clean reads) were obtained by trimming Smart-seq2 45 public primer sequences and removing reads containing adapters, reads with more than 5% unknown nucleotides, and low-quality reads from raw data. The clean data were then aligned to the reference genome using HISAT2 v2.1.0 46 with default parameters. To estimate the expression levels of genes in each sample, the read numbers for each gene in each sample were counted by HTSeq v0.6.0 47 , and then normalized to fragments per kilobase million mapped reads (FPKM). Sequencing data in this study were deposited in the NCBI Short Read Archive (SRA) under the accession PRJNA692810. Validation by quantitative real-time PCR To validate the RNA-seq data, the expression levels of 10 DEGs were evaluated by qRT-PCR using the cDNA samples used for RNA-Seq library construction as templates. Gene-speci c primers (Supplementary  Table 4) were designed using a primer-blast program from NCBI (http://www.ncbi.nlm.nih.gov/tools/primer-blast/). qRT-PCR was performed using the Roche Light Cycler 480 (Roche, Basel, Switzerland) at the following conditions: 95℃ for 30 s, followed by 45 cycles at 95°C for 10 s, 57°C for 10 s, and 72°C for 15 s. Actin was used as a housekeeping gene to normalize the ampli cation. The relative transcript levels for each sample were calculated using the 2 −∆∆Ct method.

Metabolome analyses
A total of 20 samples, including two biological replicates at ve time points, were used for metabolomic analysis. Sample preparation and UPLC-MS/MS (UPLC, Shim-pack UFLC SHIMADZU CBM30A system, www.shimadzu.com.cn/; MS, Applied Biosystems4500 Q TRAP, www.appliedbiosystems.com.cn/) was performed as described previously 48 . Quality control (QC) samples were prepared by mixing extracts from samples to analyze the repeatability. During instrumental analysis, a quality control sample was inserted into every 10 sample analyses to monitor the repeatability of the analysis process. Based on the self-built database MWDB (Metware database); the metabolites of the samples were qualitatively and quantitatively analyzed by mass spectrometry data obtained using Analyst 1.6.1 software. According to the secondary spectrum information, isotope signals, including repeated signals of K + , N + , and NH 4 + ions, and fragment ions of other larger molecular weight substances were removed during the analysis.
Signi cantly differentially accumulated metabolites (DAMs) between groups were determined by VIP ≥ 1 and absolute Log 2 FC (fold change) ≥ 1. VIP values were extracted from orthogonal projections to latent structures-discriminant analysis (OPLS-DA) results, which also contain score plots, and permutation plots were generated using the R package MetaboAnalystR 49 . The data were log-transformed (log 2 ) and meancentered before the OPLS-DA. To avoid over tting, a permutation test (200 permutations) was performed.
Bioassays and eld validation in the greenhouse The toxicity of the selected metabolites (betulin and quercetin) to GPA was measured using the slip-dip method according to Zhou et al 50 . Five hundred and forty Apterous adult GPAs were used as test insects for each compound. The compounds were dissolved in a solution of 3% acetone containing 0.1% Tween 80. Six doses (2, 1, 0.5, 0.25, 0.125, and 0.0625 g/mL) were set with three replicates. Tobacco (Nicotiana tabacum) leaf discs (1×1 cm) were added to the corresponding solutions for 5 s, removed, and dried. Tobacco leaf discs treated with a solution of 3% acetone containing Tween 80 alone were used as the blank control (CT). The aphids were observed after 24 h and 48 h of rearing under controlled growth conditions, as described above. Aphids that exhibited immobility or irregularly trembling legs were considered dead. The lethal concentrations for subsequent experiments were determined based on logprobit analysis of concentration-mortality data.
The inhibitory effect of betulin against aphids was evaluated under greenhouse conditions. S-27 shoots, sprayed with betulin (396.08 mg⋅L − 1 , 24 h LC 50 ) or acetone (as CK), were used as test plant materials.
Each treatment consisted of three shoots. Every shoot was infested with 15 adult GPAs after it was sprayed with the corresponding solutions and dried. The number of GPAs in the shoots was recorded daily for a week. The experiment was carried out at 25 ± 1℃, 60-80% RH, and a 14:10 h light/dark photoperiod.

Declarations
Data availability. All data related to this paper is available upon request.