Metabarcoding and Metagenomics : Research Article
Print
Research Article
Disentangling higher trophic level interactions in the cabbage aphid food web using high-throughput DNA sequencing
expand article infoMarie-Caroline Lefort‡,§, Stephen D. Wratten§, Antonino Cusumano|, Yann-David Varennes, Stephane Boyer#,
‡ Unitec Institute of Technology, Auckland, New Zealand
§ Bio-Protection Research Centre, Lincoln, New Zealand
| Wageningen Universitym, Wageningen, Netherlands
¶ Fondation Rurale Interjurassienne, Courtételle, Switzerland
# Insect Biology Research Institute (IRBI), UMR 7261 CNRS/Université François-Rabelais de Tours, Tours, France
Open Access

Abstract

The lack of understanding of complex food-web interactions has been a major gap in the history of biological control. In particular, a better understanding of the functioning of pest food-webs and how they vary between native and invaded geographical ranges is of prime interest for biological control research and associated integrated pest management. Technical limitations associated with the deciphering of complex food-webs can now be largely overcome by the use of high throughput DNA sequencing techniques such as Illumina MiSeq. We tested the efficiency of this next generation sequencing technology in a metabarcoding approach, to study aphid food-webs using the cabbage aphid as model. We compared the variations in structure and composition of aphid food-webs in the species’ native range (United Kingdom, UK) and in an invaded range (New Zealand, NZ). We showed that Illumina MiSeq is a well suited technology to study complex aphid food-webs from aphid mummies. We found an unexpectedly high top down pressure in the NZ cabbage aphid food-web, which coupled to a large ratio of consumer species / prey species and a lack of potential inter-specific competition between primary parasitoids, could cause the NZ food-web to be more vulnerable than the UK one. This study also reports for the first time the occurrence of a new hyperparasitoid species in NZ, as well as new associations between hyperparasitoids parasitoids and the cabbage aphid in this country. We conclude that the complexity of aphid food-webs in agricultural systems could often be underestimated, particularly at higher trophic levels; and that the use of high throughput DNA sequencing tools, could largely help to overcome this impediment.

Keywords

metabarcoding, biological control, enemy release hypothesis, hyperparasitism, parasitoids, hyperparasitoids, competition

Introduction

According to the Enemy Release Hypothesis (ERH) (Elton 1958), phytophagous insects that migrate or are introduced into a new region may be released from predation and parasitism pressure. This is because some or all of their natural enemies may not occur in the newly colonized region, either because it is outside of their distribution range or because they did not manage to migrate themselves (MacLeod et al. 2010). This scenario is supported by the island biogeography theory, which, inter alia, predicts lower species richness on the smallest islands and on those most isolated from the nearest neighbours and the mainland (MacArthur and Wilson 1967, Sax 2008). As a consequence, predation pressure may be weaker for a pest population that has invaded an island compared to conspecifics occurring in the species’ historical range. When enemies are rare or absent and food resources are abundant (as is the case in agricultural systems), major population outbreaks can result (Wallner 1987). However, if natural enemies colonized the new region as well, they may themselves be ‘freed’ from the restraints imposed by higher trophic levels after migration (e.g. Gómez-Marco et al. 2015). This would mechanistically result in higher predation pressure on the pest species and is the principle on which classical biological control hinges.

A better understanding of the functioning of these food webs and how they vary between native and invaded geographical ranges is of prime interest for biological control research and integrated pest management (Varennes 2016, Gómez-Marco et al. 2015, Gurr et al. 2017). Such knowledge is essential to design more efficient classical biological control of pests through strategic disruption of the established trophic interactions in a given geographical range. This is particularly important when pest species are involved in food webs comprising four trophic levels or more as is the case in interactions involving host plants - aphids - parasitoids and hyperparasitoids.

Aphids are widely recognised as a major pest for a wide variety of crops (Van Emden and Harrington 2007); however, the intensity of their impact varies with geographical location. For example, the cabbage aphid (Brevicoryne brassicae (L.) is a prevalent oilseed rape (Brassica napus L. (Brassicaceae)) pest in New Zealand (NZ) (Ellis and Singh 1993), and has been a problem for almost a century in this invaded region (Ellis and Farrell 1995). In contrast, in the United Kingdom (UK), which is part of its native range (Blackman and Eastop 2000), B. brassicae is only considered as locally and occasionally damaging (Ellis et al. 1999, Alford et al. 2003).

In aphid food webs, the taxonomic identification of upper trophic levels, mainly composed of minute species of parasitic wasps and flies, is often difficult because of the lack of taxonomic expertise for the genera and species of interest (Smith et al. 2011). This can lead to inaccurate identifications (Day 1994, Derocles et al. 2012). Recent progress has been made to overcome this impediment through a more widespread utilisation of molecular tools in food web studies (e.g. Valentini et al. 2009, Pompanon et al. 2011, Andrew et al. 2013, Boyer et al. 2013). Despite this recent progress, a refinement of the methodologies to study aphid-based food webs is still needed. In a recent study, Varennes et al. (2014) proposed a new molecular method for the construction of aphid-based food webs based on the amplification of aphid, parasitoid and hyperparasitoid DNA from empty mummies (i.e. after emergence) (Fig. 1). The DNA left inside aphid mummies was successfully amplified using universal primers and species could be identified by single-stranded conformation polymorphism analysis (SSCP). However, this method often failed at amplifying the DNA of multiple trophic levels at once. In the mixed sample that is an empty mummy, the DNA from one species would often be present in higher quantity and as a result preferential-amplification could occur and mask the presence of other species on the SSCP gel. In addition, the actual identification of the species involved requires a priori establishment of a library of banding patterns for all species potentially present. These limitations can be largely overcome by the use of high throughput DNA sequencing (Furlong 2014, González-Chang et al. 2016) such as Illumina sequencing, where the actual sequences of DNA are obtained for all species present in a mixed DNA samples.

Figure 1.

Opened aphid mummy (second trophic level) containing a near-adult parasitoid (third trophic level) and hyperparasitoid eggs (fourth trophic level). Modified from Varennes (2016).

We tested the efficiency of Illumina MiSeq technology to describe aphid food webs, and then compared the variations in both structure and composition of webs in the species’ native range (UK) as well as in an invaded range (NZ). In view of the relative isolation of New Zealand from the nearest mainand and its distance to the native range of the cabbage aphid, and in accordance to the island biogeography theory, we hypothesised that food webs will comprise fewer species and fewer trophic levels in the aphid’s invaded range.

Material and Methods

Insect sampling and DNA extraction

A total of 99 aphid mummies were collected from five oil seed rape fields in NZ (n=50) and in the UK (n=49) (Suppl. material 1). Mummies were collected during summer (i.e. end of November in NZ and mid-July in the UK), which corresponds to the abundance peak period of the aphid populations. The collector walked 50 m long transects, located at a minimum of three meters from the edge of the field to reduce edge effects. In the field, mummies were individually stored on 96-well microcentrifuge plates filled with propylene glycol. Plates were then shipped (UK) or transported (NZ) back to the laboratory and stored in a -80°C freezer until processing. For each of the 99 mummies, a DNA extract was prepared using the ZR Tissue and Insect DNA MicroPrep extraction kit (Zymo Research) following the protocol by Varennes et al. (2014) with the following modifications. After being crushed in a bead beater, samples were left to incubate overnight in a water bath at 55° C, and 2x20 µL of elution buffer was passed through the Zymo-Spin IC column and left for 30 min at room temperature prior to centrifugation, resulting in a 40 µL final eluted DNA solution.

DNA sequencing

Illumina MiSeq sequencing

For each sample, 1 µL of DNA extract was used as template to prepare DNA libraries prior to sequencing with Illumina MiSeq. All DNA amplifications and the subsequent Illumina MiSeq run were performed by New Zealand Genomics Limited (NZGL). DNA was amplified with fusion primers comprising universal primers designed to amplify a 455-bp internal region of the mitochondrial gene COI (i.e. MLepF1 (forward GCTTTCCCACGAATAAATAATA) and LepR1 (reverse TAAACTTCTGGATGTCCAAAAAATCA)) (Hajibabaei et al. 2006), molecular identifiers as recommended by the manufacturer, and adapters for Illumina MiSeq attached using Nextera XT Index Kit (Illumina 2011). Samples were pooled from 96 wells plates at equimolar concentration (2 ng/uL) whenever possible. Few recalcitrant samples had lower yield but were nonetheless included in the run. The 2x300bp Illumina MiSeq run comprised a total of 384 samples but only 99 of these (about 25%) corresponded to the current study.

Bio-informatic analyses

Bioinformatic analyses were performed by NZGL. The Illumina sequences were checked through the SolexaQA++ v.3.1.4 (Cox et al. 2010) analysis pipeline and sorted by country of origin. Reads that contained a base where the reported quality score dropped below 3 were trimmed, with the remainder of the sequence being discarded. Primers were trimmed and all pairs of reads were then merged using USEARCH v7 (Edgar 2010) and sequences shorter than 200 base pairs after merging were discarded. VSEARCH v2.0.3. (Rognes et al. 2016) was used to filter and cluster the resulting Illumina MiSeq data (command used: USEARCH –fastq_mergepairs –fastq_trunqual 3 –fastq_minmergelen 200), discarding all sequences that had more than one expected error per read (maxee = 1.0) (Edgar and Flyvbjerg 2015).

All data were then de-replicated (i.e. all non-unique sequences were removed, to make downstream computation faster) and all sequences that occurred only once in the overall dataset were discarded. The unique sequences were then clustered at 97% identity to form Molecular Operational Taxonomic Units (MOTUs), using the cluster_fast command in VSEARCH. Any MOTU that contained only one sequence was regarded as a potential sequencing error and was therefore discarded. This happened only once with our dataset.

Taxonomy for each cluster was assigned via a BLAST (Camacho et al. 2009) search for each unique sequence against GenBank (Benson et al. 2005) and the BOLD database. The percentage identity threshold was conservatively set at 98% to assign species name to MOTUs. MOTUs corresponding to the aphid (B. brassicae) were only used to calculate detection rates and then discarded to focus only on the analysis of parasitoids and hyperparasitoids.

Chimeric sequences, PCR artefacts and missing information on the databases could all lead to low percentage identity, therefore, for MOTUs with percentage identity below the 98% threshold, only those that stricly produced more than five reads in at least one sample and were detected in at least two samples were considered robust and retained in the food web analysis (Boyer et al. 2012, Waterhouse et al. 2014). The retained MOTUs are assumed to represent species that are absent from the Genbank and BOLD databases. In this case, the closest match was used to determine whether the MOTU was a parasitoid or a hyperparasitoid. Because of the degraded state of the environmental DNA used here, it is possible that DNA alterations would have generated 'erroneous' sequences potentially leading to unidentified MOTU. Therefore, unidentified MOTU should be treated with caution with regards to their reality in nature.

Trophic levels detection and food webs composition

The detection rates of the different trophic levels in each food web were compared using Mc Nemar statistical tests (α = 0.05).

The species composition of each trophic level was determined via BLASTn searches (Karlin and Altschul 1990, Camacho et al. 2009) on the GenBank and BOLD databases as described in the Bio-informatic analyses section. The structure of the two food webs was described and compared using standard food web connectedness descriptors, namely taxa richness, trophic links, connectance, ratio of prey to consumer and food web vulnerability (Thompson et al. 2012). Hyperparasitism rates (i.e. ratio of hyperparasitoid to mummy) in the two study sites (i.e. NZ and the UK) were compared with a Chi square test, using the statistical software R (R Core Team 2014).

Results

Detection of trophic levels

All samples were successfully sequenced using Illumina MiSeq technology. A total of 672,941 merged reads were obtained from NZ samples and 803,676 from UK samples. When combining reads from both countries, the detection rate for the third trophic level (i.e. parasitoid) reached 99%, while only just over 37% of the samples analysed with Illumina MiSeq produced DNA sequences for the second trophic level (i.e. aphid) (Fig. 2), which constitutes a highly significant difference (Mc Nemar test, X2 = 8.0497, df = 1, P < 0.01). DNA from the fourth trophic level (i.e. hyperparasitoid) was detected in almost 74% of the samples.

Figure 2.

Comparison of DNA detection rate of three trophic levels in Brevicoryne brassicae food web using Illumina MiSeq technology. ** indicates P<0.01 for Mc Nemar test at 5% significance level.

Third and fourth trophic level composition

The UK food web comprised three parasitoids: two identified species, Diaeretiella rapae (McIntosh) and Aphidius rhopalosiphi de Stefani-Perez and one unidentified MOTU. There were also two species of hyperparasitoid: Alloxysta leunisii (Hartig) and Al. tscheki (Giraud) (Fig. 3). The New Zealand food web comprised only one parasitoid: D. rapae, but eight potential hyperparasitoids: Alloxysta consobrina (Zetterstedt), Al. tscheki, Asaphes vulgaris Walker and four unidentified MOTUs (Fig. 3). All sequences, MOTU identification, assigned taxonomy, and occurrence (i.e. presence / absence) in individual mummies as well as exploratory statistics addressing sequencing depth per country and MOTU rarefaction are available as supplementary files (Suppl. materials 2, 3).

Figure 3.

Brevicoryne brassicae food web structure in native (UK) and invaded (NZ) ranges, and connectedness descriptors (see Bersier et al. 2002 and Thompson et al. 2012), based on the molecular study of aphid mummies (N=99). Full circles represent trophic levels identified to species level based on COI DNA sequences (>98% BLAST similarity) and full circles with disrupted borders represent MOTUs which could not be identified down to species based on the Genbank and BOLD databases.

Food web complexity, direct connectance and prey/consumer ratio were higher in the UK food web, while vulnerability was higher in the New Zealand food web (Fig. 3).

The overall hyperparasitism rates did not significantly differ between the native (UK) and the invaded (NZ) ranges (Chi2 = 0.5759, df = 1, p = 0.4479). However a much higher species richness was apparent at the fourth trophic level in the New Zealand food web (Fig. 4).

Figure 4.

Interacting species within individual mummies of the aphid Brevicoryne brassicae collected from native range (UK) and invaded range (NZ), depicting inter-specific interactions in the upper tropic levels.

Inter-specific competition within the fourth trophic level was detected in both food webs in the form of multiple species of hyperparasitoid present in an individual mummy (Fig. 4). There was also a strong inter-specific competition at the third trophic level in the UK food web, where the DNA of more than one parasitoid was retrieved in over 75% of the mummies. However, the competition appears higher at the fourth trophic level in the NZ food web, where 30 % of the mummies contained the DNA of more than one hyperparasitoid. Only two species were sometimes found as sole hyperparasitoid (namely Al. consobrina in NZ and Al. leussini in the UK). All other hyperparasitoid taxa were always found in the presence of one of these two species.

Discussion

Amplification success and detection of trophic levels

In the present study, the use of Illumina MiSeq technology allowed the detection of at least three trophic levels from the food web studied (i.e. trophic levels 2-4). Although older technologies such as 454 pyrosequening presented the advantage of producing longer targeted amplicons (Dowle et al. 2015), this asset is not a requirement for aphid food web reconstruction, where short sequences of degraded DNA are usually targeted (Varennes et al. 2014) as they are sufficient to complete species detection (Meusnier et al. 2008).

In aphid food webs, it can be assumed that as the parasitoid larva develops and consumes its host, the amount of aphid DNA inside the mummy decreases while that of parasitoid DNA increases. After parasitoid emergence, the quantities of leftover DNA from the host itself are therefore likely to be low compared to that of parasitoids, which could explain the rather low detection success of aphid DNA from mummies (i.e. less than 38%). The detection of aphid DNA may have also been affected by primer bias (Piñol et al. 2014, Elbrecht and Leese 2015) potentially leading to preferential amplification of DNA from higher trophic levels. Because the aim of the present study was primarily descriptive (i.e. food web reconstruction), potential quantification bias had little bearings on the results. Besides, DNA quality is expected to quickly degrade in the field through the action of microbial metabolism and exonuclease activity (Strickler et al. 2015), particularly at the time of sampling (i.e. summer) when temperatures are at their highest. In this study, mummies were collected from fields regardless of their age, which might further explain the general low detection rates of aphid DNA. However, from a practical point of view, this should not be regarded as a real impediment for the study of aphid food webs since mummies will either be produced in the laboratory from known species of aphids or collected from the field from aphid colonies easily identifiable via traditional taxonomic or molecular methods.

In contrary to the second trophic level, the detection rate was quite high for the third trophic level, with every single sample producing parasitoid sequences. It has been previously shown that parasitoid and hyperparasitoid DNA can be easily retrieved from aphid mummies in which the parasitoid or hyperparasitoid has developed (Traugott et al. 2008, Derocles et al. 2012, Gariepy et al. 2013). In empty mummies, a large part of parasitoids (and hyperparasitoids) residual DNA may come from their faecal pellet, meconium (Haviland 1922, El-Heneidy and Adly 2015) and/or exuviae (Lefort et al. 2012). Although not recorded in this study, it appeared that the quantities of parasitoid and/hyperparasitoid excrement can greatly vary from one mummy to another (personal observation). Selecting mummies that contain large amounts of parasitoids/hyperparasitoids excrements in aphid food web studies might significantly help to improve successful DNA amplification rate.

Hyperparasitism is considered to be one of the major causes of mortality in aphid primary parasitoids (Chua 1977, Rosenheim 1998, Sullivan and Völkl 1999, Gómez-Marco et al. 2015). To our knowledge, the highest hyperparasitism rate of aphids in oilseed rape cultures was reported to be around 50 % during the late growing season (Varennes 2016), and generally appears lower than this the rest of the time on this crop (e.g. in Iran, see Nematollahi et al. 2014). Regarding the aphid B. brassicae higher rates of up to 76.1% have also been reported elsewhere (Bahana and Karuhize 1986). In the present study, it was not possible to confidently evaluate the detection success of hyperparasitoids since there was no mean to externally determine whether the mummies used had been hyperparasitised, and samples were not dissected prior to destruction for sequencing purpose to avoid contamination. There are several types of aphid hyperparasitoids, those which attack the aphid before mummification is completed and lay their eggs inside the body of the primary parasitoid larva, those which attack the mummy “irrespective of whether it contains primary or a secondary parasitoids” (Müller et al. 1999), and species which use a dual oviposition strategy (i.e. attacking either living aphids or mummies) such as the aphid hyperparasitoid Syrphophagus aphidivours (Sullivan and Völkl 1999). In either case, considering that parasitoids detection level reached 100% with Illumina, despite parasitoid residual DNA likely to have been older than that of hyperparasitoids, it is assumed that all (or most) hyperparasitoid DNA would have also been amplified successfully as their DNA is likely to be fresher and better preserved than that of parasitoids. While bearing in mind that this method is based on DNA detection only, and do not distinguish between oviposition success and successful development to a mature hyperparasitoid, the DNA amplification rate for hyperparasitoids might hence be a good indication of the true hyperparasitism rates occurring in the studied populations of cabbage aphid. With 63.3% in the UK food web and 84% in the NZ food web.

The results of the present study demonstrated that, with 100% of DNA amplification success and high detection rates of DNA from multiple trophic levels, Illumina MiSeq is a suitable technology to study aphid food webs composed of more than two trophic levels.

Food webs’ structure

The molecular study of B. brassicae food webs, revealed a much more complex structure in the aphid’s invaded range (i.e. NZ) compared to its native geographical range (i.e. UK). The NZ food web appeared to contain fewer primary parasitoids, which is concordant with the ERH (Elton 1958) and the island biogeography theory (MacArthur and Wilson 1967, Sax 2008). However, the NZ food web also contained an increased number of hyperparasitoids species, which is not consistent with the island biogeography theory, as lower species richness would be expected (MacArthur and Wilson 1967, Sax 2008).

In the current study, samples were collected during what was considered the peak abundance period, and while the UK hyperparasitism rate already appears quite high (i.e. 63.3%), the equivalent rate for the NZ mummies (i.e. 84%) dwarfs all previous reports (see Bahana and Karuhize 1986, Nematollahi et al. 2014, Varennes 2016). While, these rates did not statistically differ between the two countries, the top-down pressure could be exacerbated in the invaded range due to higher species richness in hyperparasitoids. While the pressure imposed by one or two hyperparasitoids might be limited in time (Nematollahi et al. 2014), the pressure imposed by eight different species (as observed in NZ) might overlap over longer periods of time, and significantly weaken biological control by parasitoids. Indeed, intra-seasonal fluctuations by hyperparasitoid genera and/or species have been reported in many aphid food web studies. For instance, in a study published by Lohaus et al. (2012), the abundance and composition of the fourth trophic level varies across the season, where two hyperparasitoid genera, Alloxysta and Phaenoglyphis, were essentially represented during the flowing period and where three other genera increased density during the peak ripening of conventional wheat fields. Similarly, Gómez-Marco et al. (2015) observed an increased abundance of hyperparasitoids belonging to the genus Alloxysta at the beginning of the season versus an increased abundance of Syrphophagus aphidivorus (Mayr) (Encyrtidae) toward the end of the season in the food web of Aphis spiraecola Patch in clementine orchards. In 2012, Gurr et al. (2012) stressed that the main reason for the low success rate of biological control in NZ was the failure of the biological control agent to establish. The suggested increased top-down pressure imposed by multiple species at the upper trophic levels could also explain why the success of biological control in NZ tends to remain so low.

In addition to the increased top down pressure described above, the large ratio of consumer species (parasitoids + hyperparasitoids) / prey species (parasitoids + aphid) and the lack of potential for inter-specific competition between primary parasitoids, causes the NZ food web to be more vulnerable than the UK one (i.e. VNZ = 3.33, VUK = 2). As a result, the third and fourth trophic levels of the system in the invaded range (i.e. NZ) are more likely to collapse and the associated biological control to fail. Furthermore, the third trophic level of the NZ food web appeared to be only composed of D. rapae, which renders this web even more vulnerable and subject to collapse. It is important to note that the conservative 2% species delineation threshold chosen for this study limits unnecessary incorporation of wrongly identified species in the food web. However, other food web studies sometimes report the use of higher thresholds (e.g. Prévot et al. 2013). Puillandre et al. (2011) insisted on the importance of making species delimitation more reliable in DNA barcoding studies, and care must be taken with interpretation based on these thresholds.

Food webs’ composition

Varennes (2016) suggested potential high differences in species composition in food web of OSR between New Zealand and Europe, where this plant has been cultivated for more than five centuries (Bunting 1985). The present study, seems to confirm this hypothesis since eight different species of hyperparasitoids were detected in the NZ food web, among which only one also occurs in the UK food web for this aphid species. Our findings also confirm that molecular tools have the potential to unravel higher species richness at fourth trophic level; food web studies usually only report two or three species in aphid webs, regardless of the host plant (e.g. Lotfalizadeh 2002, Vaz et al. 2004, Nematollahi et al. 2014). Nevertheless, because not all MOTUs could be identified at species level, it is difficult to determine to what extent the food web has recruited local species in the invaded range. This impediment certainly holds on a lack of species record on the Genbank and BOLD databases. There are five Alloxysta wasps currently known to occur in NZ (Ferrer-Suay et al. 2012), and only two of them are known to be secondary parasitoids on B. brassicae (see Ferrer-Suay et al. 2014). Only three introduced hyperparasitoid species could be named with certainty (As. vulgaris, Al. consobrina and Al. tscheki, percentage similarity >98%). Other hyperparasitoid MOTUs could not be named and may be native or introduced hyperparasitoids for which no COI sequences have yet been deposited in Genbank or BOLD. Another MOTU that could not be nammed is MOTU 39 a parasitoid in the UK web. This MOTU is very closely related to D. rapae but in one region of its COI sequence there are 12 Gs, while in the reference sequence of D. rapae all these 12 bases are Ts (MOTU39: GGGGGAGGAGGAAGGCCAGCGGG, D. rapae: TTTTGATTATTAATTCCATCTTT). The explanation for this extraordinarry pattern is unknown, but the fact that MOTU 39 was the fifth most common MOTU in the UK and was found in 35 different samples (out of 49), makes it unlikely to be a sequencing artefact.

A number of native hyperparasitoid species have recently been described from NZ including the first endemic charipines for New Zealand: Al. rubidus n. sp. and Al. thorpei n. sp (Ferrer-Suay et al. 2012). However, it is not yet known whether these two species are present in agricultural landscapes (Varennes 2016). It is also important to note that our knowledge of global parasitoid diversity sensu lato might be as low as 1% (La Salle and Gauld 1992, Smith et al. 2011), which means that many species remain to be described. This is consistent with the high number of MOTUs with no species match detected at the fourth trophic level in the NZ food web.

The molecular-based identification also revealed the presence of Al. tscheki, an Asiatic charipine wasp that had never been intercepted nor recorded in New Zealand (Dave Voice, Senior Scientist, Ministry for Primary Industries, pers. com.). Our results also revealed new associations as the hyperparasitoid Al. leunissii (found in more than 60% of NZ mummies) had never been reported before as attacking D. rapae or any other parasitoid within a B. brassicae host, (see complete list of known associations in Ferrer-Suay et al. 2014). This new association could contribute to preventing D. rapae from reaching high population densities in NZ, which may counteract the full pest-suppressive potential of this biological control agent. The risk is particularly accrued by the fact that D. rapae is often reported as the dominant primary parasitoid of B. brassicae (e.g. Pike et al. 1999) and even sometimes the only one (see Nematollahi et al. 2014), particularly in brassica crops (Pike et al. 1999).

Trophic levels

A number of samples (aphid mummies) contained DNA from more than one parasitoid and/or more than one hyperparasitoid. This could be explained by competition, where several species of parasitoid or hyperparasitoids attack the same host. Such intra-guild competition has commonly been reported in parasitic wasps (Cusumano et al. 2016). In the case of hyperparasitoids, competition may be particularly important in UK because we detected only Alloxysta species which are closely related and use the same developmental strategy. They are all koinobiont endophagous, i.e. they oviposit inside a parasitoid larvae when the aphid is still alive (Sullivan and Völkl 1999). The importance of competition among natural enemies in shaping species coexistence and community structure has often been investigated from a theoretical perspective (May and Hassell 1981, Comins and Hassell 1996, Bonsall and Hassell 1999), although there have been disagreements about its importance under natural conditions (Force 1974, Dean and Ricklefs 1979, Force 1980). The potential of molecular tools to disentangle foodweb complexity can significantly advance our understanding of the role played by interspecific competition in basic and applied ecology.

Another potential explanation to the presence of DNA from multiple hyparasitoids in the same mummy, could be higher trophic relationships. While the interactions between the MOTUs were conservatively represented in three different trophic levels in this study, it is important to bear in mind that certain hyperparasitoid species can also attack other hyperparasitoid within mummies (Matejko and Sullivan 1984). Such strategy has been reported for Asaphes and Dendrocerus species which are idiobiont ectophagous, i.e. they oviposit on parasitoid prepupae and pupae inside the mummy. For example, Sullivan (1972) reported that As. californicus (Provancher), successfully attacked, oviposited on and emerged from the other hyperparasitoid Al. victrix (Westwood). By engaging in such interspecific tertiary parasitism, these species may extend the food web to a 5th trophic level. The concurrent presence of species of both Alloxysta and Asaphes genera in the NZ B. brassicae food web, and the prevalence of multiple hyperparasitoid species per mummies (54% of samples), strongly support the possibility of the existence of this 5th trophic level in this country.

In a recent study, where multiplex PCR was used to describe the food web of Aphis spiraecola Patch, Ferrer-Suay et al. (2012) reported the commonness of multiple hyperparasitoid species occurrence in a single mummy. Our study also illustrates the utility of molecular tools to describe aphid food webs, since the detection of multiple competing species at higher trophic level, or that of a 5th trophic level, could not be achieved using traditional taxonomy methods. This is because, while using the emergence method, only the one species that will successfully achieve development and emerge will be identified. However, the disentangling of which hyperparasitoid species are competing and which ones may be specifically targeting other hyperparasitoid species remains difficult. One indicator may be the systematic association of a particular species with another one. In our study, the majority of hyperparasitoids detected were never found alone. In fact, only Al. consobrina in NZ and Al. leussini in the UK, could be found as sole hyperparasitoid in a mummy. All other hyperparasitoid species were always found in association with the two aforementioned taxa. This may indicate the existence of fifth tropic level species targeting Al. consobrina in NZ and Al. leussini in the UK.

Conclusion

The complexity of food webs in agricultural systems has often been underestimated (González-Chang et al. 2016) and the present study of the B. brassicae food web seems to confirm this trend. Based on the results of the present study, we recommend the use of high throughput DNA sequencing tools, particularly Illumina sequencing, to study aphid food webs in agricultural ecosystems. Unlike multiplex PCRs, which are often used in aphid food web studies (e.g. Traugott et al. 2008, Gómez-Marco et al. 2015) and the success of which relies on a thorough knowledge of the interacting species, high throughput DNA sequencing tools also allow the discovery of new species and new species associations. In a context of climate change which has the potential to highly disrupt the structure and the functioning of local food webs and the outbreak of pests, it is essential that the tools used to understand these food webs allow us, as much as possible, to get a complete picture of the scene. This holds particularly true in invaded ranges such as NZ, where, as for B. brassicae, the food web vulnerability is increased.

Data accessibility

Raw sequence data is available on the NCBI repository under the following accession numbers, UK sequences: SRR6039612 - SRR6039660; NZ sequences: SRR6039544 - SRR6039593

Summary tables containing the number of reads for each MOTU in each sample as well as R code for data exploration are available on Figshare (DOI: 10.6084/m9.figshare.5350729).

Acknowledgements

The authors would like to thank the Bio-Protection Research Centre for its financial support. The authors also would like to thank Annie Barnes for her technical assistance in the field, as well as Richard Harrington, Matt Skellern and Sam Cook for their help with the UK samples.

Author contributions

MCL - performed the molecular analysis, analysed the data, wrote the manuscript, prepared the figures, prepared revisions of the drafts, performed the revisions

SB - performed the molecular analysis, analysed the data, reviewed drafts of the manuscript, prepared revisions of the drafts, performed the revisions

AC - reviewed drafts of the manuscript

YDV - collected the samples, reviewed draft of the manuscript

SDW - provided financial contribution to reagents/materials, reviewed draft of the manuscript

References

Supplementary materials

Suppl. material 1: Supporting Information 1 
Authors:  Lefort M-.C., Wratten S.D., Cusumano A., Varennes Y-.D., Boyer S.
Data type:  Sampling / Amplification success
Brief description: 

OSR aphid mummy collection. Sampling location and size / Amplification success of mummies’ DNA extracts by Illumina sequencing.

Suppl. material 2: Supporting Information 2 
Authors:  Lefort M.-C.
Data type:  COI sequences, MOTU identification, assigned taxonomy, and occurrence (i.e. presence / absence) in individual mummies
Suppl. material 3: Supporting information 3 
Authors:  Boyer S.
Data type:  Graphs exploring the data
Brief description: 

Exploratory statistics addressing sequencing depth per country and MOTU rarefaction.