Rapid transcriptional plasticity of duplicated gene clusters enables a clonally reproducing aphid to colonise diverse plant species

Abstract

Background: The prevailing paradigm of host-parasite evolution is that arms races lead to increasing specialisation via genetic adaptation. Insect herbivores are no exception and the majority have evolved to colonise a small number of closely related host species. Remarkably, the green peach aphid, Myzus persicae, colonises plant species across 40 families and single M. persicae clonal lineages can colonise distantly related plants. This remarkable ability makes M. persicae a highly destructive pest of many important crop species.

Results: To investigate the exceptional phenotypic plasticity of M. persicae, we sequenced the M. persicae genome and assessed how one clonal lineage responds to host plant species of different families. We show that genetically identical individuals are able to colonise distantly related host species through the differential regulation of genes belonging to aphid-expanded gene families. Multigene clusters collectively upregulate in single aphids within two days upon host switch. Furthermore, we demonstrate the functional significance of this rapid transcriptional change using RNA interference (RNAi)-mediated knock-down of genes belonging to the cathepsin B gene family. Knock-down of cathepsin B genes reduced aphid fitness, but only on the host that induced upregulation of these genes.

Conclusions: Previous research has focused on the role of genetic adaptation of parasites to their hosts. Here we show that the generalist aphid pest M. persicae is able to colonise diverse host plant species in the absence of genetic specialisation. This is achieved through rapid transcriptional plasticity of genes that have duplicated during aphid evolution.

Keywords: Gene duplication; Genome sequence; Hemiptera; Myzus persicae; Parasite; Plasticity; RNA interference (RNAi); Sap-feeding insects; Transcriptome.

Fig. 1

High rate of lineage-specific gene accumulation in aphids relative to all other insect orders. Figures show arthropod phylogenetic relationships, per genome proportions of single copy (blue) and duplicated (red) genes and orthology relationships among arthropod genes based on gene family clustering with MCL [30]. Phylogenetic relationships among arthropod species included for gene family clustering were estimated using RAxML [31] based on a protein alignment of 66 single-copy orthologs found in all taxa. This topology and protein alignment was then used to infer relative divergence times with RelTime [32] under an LG substitution model. Inset shows relative rate of lineage-specific gene accumulation for all included insect orders and comparison with aphids. Error bars show standard deviation of species within a given grouping. Relative rates of lineage-specific gene accumulation were calculated for each species by dividing the number of group specific genes (either order-specific or aphid-specific) by the crown plus stem age for the given group (in relative divergence time)

Fig. 2

M. persicae experienced greater gene loss rates (a) and stronger purifying selection in retained ancestral duplicates (b) than A. pisum. a Age distribution of duplicated genes in M. persicae and A. pisum. The number of synonymous substitutions per synonymous site (d _S) was calculated between paralog pairs for M. persicae (green) and A. pisum (blue) using the YN00 [91] model in PAML [82]. For each duplicated gene, only the most recent paralog was compared. Pairwise d _S was also calculated for 1:1 orthologs between M. persicae and A. pisum (red), the peak in which corresponds to the time of speciation between the two aphid species. After filtering, 1955 M. persicae paralog pairs, 7253 A. pisum paralog pairs and 2123 1:1 orthologs were included for comparison. Mean d _S of 1:1 orthologs between A. pisum and M. persicae was 0.26. b Box plots showing median d _N /d _S for A. pisum and M. persicae paralog pairs that duplicated before and after speciation of the two aphid species and for 1:1 orthologs between the two species. Older duplicate genes have lower d _N /d _S than recently duplicated genes (since speciation) indicating stronger purifying selection in ancestral versus recent duplicates. Additionally, older duplicate genes in M. persicae have significantly lower d _N /d _S than in A. pisum (Mann–Whitney U = 1816258, M. persicae: 1348 paralog pairs, A. pisum: 3286 paralog pairs, p = < 0.00001) indicating stronger genome streamlining in M. persicae than in A. pisum. Box plots are shaded by species as in (a)

Fig. 3

The set of differentially expressed genes of M. persicae clone O reared on B. rapa and N. benthamiana is enriched for (a) genes belonging to gene families with known functions, (b) tandemly duplicated genes in the M. persicae genome, (c) genes belonging to gene families expanded in aphids or unique to aphids, (d) duplicated genes before M. persicae and A. pisum diverged and (e) genes with stronger purifying selection than the genome-wide average. a–c Volcano plots of differentially expressed genes of M. persicae reared on B. rapa and N. benthamiana. Negative log₂ fold changes indicate upregulation on B. rapa and positive values indicate upregulation on N. benthamiana. a Differentially expressed genes from four gene families that have the highest number of differentially expressed genes are highlighted. These are: RR-2 cuticular proteins (n = 22), cathepsin B (n = 10), UDP-glucosyltransferase (n = 8) and cytochrome P450 (n = 5). b The set of differentially expressed genes is enriched for tandemly duplicated genes. c The set of differentially expressed genes is enriched for genes from families that are either significantly expanded in aphids compared to other arthropods (binomial test, main text) or are unique to aphids. d Time since most recent duplication (measured as d _S) for all paralogs in the M. persicae genome compared to those differentially expressed upon host transfer. Duplicated genes implicated in host adjustment (at least one of the pair differentially expressed) have a significantly different distribution to the genome wide average (p < 0.05, permutation test of equality) and are enriched for genes that duplicated before M. persicae and A. pisum diverged. e d _N/d _S distribution for duplicated genes differentially expressed upon host transfer vs. the genome wide average. Duplicated genes involved in host adjustment are under significantly stronger purifying selection than the genome wide average (median d _N/d _S = 0.2618 vs. 0.3338, Mann–Whitney U = 105,470, p = 1.47 × 10⁻⁴, two-tailed)

Fig. 4

Cathepsin B genes that are differentially expressed upon M. persicae host change belong predominantly to a single aphid-expanded clade and form gene clusters in the M. persicae genome. a Maximum likelihood phylogenic tree of arthropod cathepsin B protein sequences. The sequences were aligned with Muscle [76] and the phylogeny estimated using FastTree [92] (JTT + CAT rate variation). Circles on branches indicate SH-like local support values >80%, scale bar below indicates 0.1 substitutions per site. Rings from outside to inside: ring 1, M. persicae cathepsin B (MpCathB) gene identities (IDs) with numbers in red indicating upregulation of these genes in M. persicae reared for one year on B. rapa relative to those reared for one year on N. benthamiana and bold font indicating location on the cathepsin B multigene clusters shown in (b); ring 2, red squares indicating MpCathB genes that are differentially expressed upon M. persicae host change; ring 3, cathB genes from different arthropods following the colour scheme of the legend in the upper left corner and matching the colours of the branches of the phylogenetic tree; ring 4, aphid-expanded (AE) clades with AE_Clade I labelled light green and AE_Clade II light blue. b MpCathB multigene clusters of the M. persicae genome. Lines indicate the genomic scaffolds on which the MpCathB genes are indicated with block arrows. Gene IDs above the genes match those of the phylogenetic tree in A, with block arrows and fonts highlighted in red being differentially expressed upon host change. Scale bar on right shows 1 kb. c Relative expression levels of MpCathB genes of M. persicae at seven weeks being reared on N. benthamiana (Nb), B. rapa (Br) and A. thaliana (At). Numbers under the graphs indicate MpCathB gene IDs with those in red font DE as in (a). Batches of five adult females were harvested for RNA extraction and quantitative real-time polymerase chain reaction assays. Bars represent expression values (mean ± standard deviation (SD)) of three independent biological replicates. *p < 0.05 (ANOVA with Fishers LSD to control for multiple tests). d As in (c), except that individual aphids reared on At were transferred to At (At to At) or Nb (At to Nb) and harvested at two days upon transfer. e As in (d), except that individual aphids reared on Nb were transferred to Nb (Nb to Nb) or At (Nb to At) and harvested at two days upon transfer

Fig. 5

RNAi-mediated knock-down of the expression of multiple cathepsin B genes reduces M. persicae survival and fecundity on A. thaliana. a Relative cathepsin B (CathB) expression levels (compared to aphids on dsGFP (control) plants) of M. persicae on three independent transgenic lines (lines 5–1, 17–5 and 18–2) producing double-stranded (ds) RNA corresponding to multiple M. persicae cathepsin B genes (dsCathB) (Fig. 3a, Additional file 20: Figure S10). Aphids were reared on the transgenic lines for four generations. Batches of five adult females were harvested for RNA extraction and qRT-PCR assays. Bars represent expression values (mean ± standard deviation (SD)) of three independent biological replicates. b CathB-RNAi M. persicae produces less progeny compared to control (dsGFP-treated) aphids on A. thaliana. Five nymphs were transferred to single plants and produced nymphs on approximately day 5. Nymph counts were conducted on days 7, 9 and 11 and removed. Columns show the mean ± SD of the total nymph counts for these three days of three biological replicates, with each replicate consisting nymphs produced by 15 aphids at five aphids per plant (n = 3 plants). c, d Survival rates of CathB-RNAi and control (dsGFP-exposed) M persicae on non-transgenic A. thaliana (At) and N. benthamiana (Nb) plants. Ten third instar nymphs on dsCathB and dsGFP transgenic plants were transferred to non-transgenic plants; survival rates were recorded two days later. Bars represent mean ± SD of three biological replicates, with each replicate consisting of the survival rates of 30 aphids at 10 aphids per plants (n = 3 plants). e, f Fecundity rates of CathB-RNAi and control (dsGFP-exposed) M. persicae on non-transgenic A. thaliana (At) and N. benthamiana (Nb) plants. Nymph counts were conducted as in (b). Asterisks (*) and different letters (a, b) above the bars indicate significant difference at p < 0.05 (ANOVA with Fisher’s LSD to control for multiple tests)