Continental-scale suppression of an invasive pest by a host-specific parasitoid underlines both environmental and economic benefits of arthropod biological control

Multi-country pest & natural enemy survey

From early 2014 until late 2017, insect surveys were carried out in 601 cassava fields in Myanmar, Thailand, Lao PDR, Cambodia, Vietnam, southern China and Indonesia. Survey protocols are described in detail in Graziosi et al. (2016). In brief, we selected older fields (i.e., 8–10 months of age) in the main cassava-growing areas of each country, with individual sites located at least 1 km apart. Five linear transects were randomly chosen per site, with ten plants (routinely spaced at 0.8–1.2 m) sampled in each transect. By doing so, a total of 50 plants per field were assessed for P. manihoti infestation and per-plant mealybug abundance. In-field identification of mealybugs was based on morphological characteristics such as coloration and presence or length of abdominal waxy filaments, while samples were also taken to the laboratory for identification by specialist taxonomists. Following transect walks we calculated average P. manihoti abundance (i.e., number of individuals per infested tip) and field-level incidence (i.e., proportion of P. manihoti-infested tips per field).

To assess local A. lopezi establishment and parasitism rates, we conducted dry-season sampling from 2014 to 2017 at sub-sets of mealybug-invaded sites in Thailand (n = 20), Cambodia (n = 10, 15), southern Vietnam (n = 20, 20, 6) and Indonesia (n = 10, 9, 21) (total n = 131). Sampling consisted of collecting 20 mealybug-infested tips from local fields and transferring them to a laboratory to monitor subsequent parasitoid emergence (Neuenschwander et al., 1989). Surveys were carried out during January-May 2014 (dry season), October–November 2014 (late rainy season), January–March 2015 (dry season) in mainland Southeast Asia, and during October–November 2014 and 2017 (dry season) in insular Indonesia. Locations were recorded using a handheld GPS unit (Garmin Ltd, Olathe, KS). In-field identification of mealybugs was based on morphological characters, while samples were also transferred to the laboratory for further taxonomic identification. Voucher specimens of P. manihoti were equally deposited at the Thai Department of Agriculture (Bangkok, Thailand), Bogor Agricultural University (Bogor, Indonesia) and Plant Protection Research Institute (Hanoi, Vietnam).

To assess local A. lopezi establishment and parasitism rates, mealybug-infected tips were collected in the field and transferred to a laboratory. Upon arrival in the laboratory, each tip was carefully examined, predators were removed and the total number of P. manihoti was determined. Tips were then placed singly into transparent polyvinyl chloride (PVC) containers, closed with fine cotton fabric mesh. Over the course of three weeks, containers were inspected on a daily basis for emergence of parasitoids and A. lopezi parasitism levels (per tip and field) were computed. Next, for fields where presence of A. lopezi was reported, we carried out a regression analysis to relate field-level mealybug abundance with parasitism rate. Mealybug infestation levels and parasitism rates were log-transformed to meet assumptions of normality and homoscedasticity, and all statistical analyses were conducted using SPSS.

Multi-year mealybug and parasitoid population assessment in Vietnam

From July 2013 until July 2015, we conducted population surveys in Tay Ninh province, Vietnam; an area with near-continuous, all-year cassava cultivation (see also Le et al., 2018). The cassava mealybug is assumed to have arrived in southern Vietnam during 2011–2012, and A. lopezi was first detected from Tay Ninh province in early 2013. Eight newly-planted cassava fields were selected of uniform age, crop variety, developmental stage and management. Every two months, insect surveys were done within these fields to characterize P. manihoti incidence, infestation pressure and A. lopezi parasitism rate. In each field, a total of five linear 10–15 m transects were screened (plants routinely spaced at 0.8–1.2 m) and, 50 plants were carefully inspected for P. manihoti. Phenacoccus manihoti infestation was recorded as field-level abundance (number of individuals per infected tip) and field-level incidence (proportion of mealybug-affected tips) at each sampling date and location. To assess A. lopezi parasitism rates, 20 mealybug-infested tips were randomly collected from each field by breaking off the top parts of individual plants, and transferred to the laboratory. Parasitism rates were estimated from these samples as described above, and parasitism levels were computed for each individual field and sampling date. We used analysis of variance (PROC MIXED, SAS version 9.1; SAS Institute, Cary, NC, USA) with field as the random factor, and tested the effect of cassava age, sampling date and year for P. manihoti incidence, abundance and A. lopezi parasitism. Means were compared with least squares means approach. Mealybug abundance data were log-transformed while incidence, parasitism and hyperparasitism data were arcsine-transformed to meet normality.

The intrinsic rate of mealybug population increase, r, over 2 months was calculated over subsequent sampling events as ln(m_t+1/ m_t) where m = the per-tip mealybug density. This growth rate was regressed against the mealybug parasitism rate as a means of evaluating the role of the parasitoids in suppressing mealybug population growth rates and also of estimating the parasitism level that correlates with a decrease in mealybug population growth rate (for other uses of this approach, see Lin & Ives, 2003; Plecas et al., 2014). The statistical significance of the relationship between parasitism rate and mealybug population growth was assessed using a generalized linear model incorporating normal error distribution with r as the response variable and parasitism level and field identity as independent variables.

Exclusion cage assays

In August 2014, a field study was initiated at the Rayong Field Crops Research Center (RYFCRC) in Rayong, Thailand (see Thancharoen et al., 2018). To assess the relative contribution of natural enemies such as A. lopezi to pest control, we employed exclusion assays (Snyder & Wise, 2001; Costamagna, Landis & DiFonzo, 2007). More specifically, to determine separate and joint effects of P. manihoti and A. lopezi on cassava crop yield, four different treatments were established using two common cassava clones: Kasetsart 50 (KU50) and Rayong 72 (R72). Treatments consisted of the following: (1) ‘full cage’ assays, in which a plant was entirely covered by a mesh screen cage to exclude all natural enemies; (2) ‘sham’ cage assays, in which a plant was covered by a screen cage to provide a microhabitat similar to that of the ‘full cage’, but left open at the sides to allow natural enemy access; (3) ‘no cage’ assays, in which a plant was kept without a cage as a ‘real-world’ benchmark. For the ‘no cage’ assay, two different mealybug infestation modes were adopted: (1) one single infestation with 10 mealybug adults per plot at the onset of the experiment; (2) a monthly infestation with 10 P. manihoti adults per plot, to ensure a sustained mealybug population over the entire course of the experiment. Each treatment (i.e., a single cage containing four replicated plants) was established with four replicates, for two different cassava clones. The experimental field was established using locally-sourced stem cuttings of KU50 or R72, planted at 1-m distances within plots.

Once plants had reached 4.5 months of age, 2 × 2 × 2 m polyvinylchloride (PVC) frame cages were deployed, with four plants contained within each cage. Cages were covered with fine nylon mesh screen to prevent entry by insects, including A. lopezi parasitoids. In January 2015, 10 adult female P. manihoti were gently brushed onto plants within each treatment. Mealybug adults were obtained from a laboratory colony at RYFCRC that had been started in early 2014, in which P. manihoti were maintained on potted cassava plants within a screen-house that were regularly supplemented with field-collected individuals. Visual observations were carried out within the cages on a monthly basis and P. manihoti abundance was recorded on each plant. On September 7, 2015, once the crop had reached 12 months of age, cages were removed and plants within the different experimental treatments were harvested manually. At harvest, fresh root yield (FRY) was determined for each plant to determine treatment effects (Karlström et al., 2016).

Mealybug population build-up under each experimental treatment was calculated, by converting the average number of mealybugs per plant (averaged across plants within a given treatment, as to avoid pseudo-replication) on a given sampling date to cumulative mealybug-days (CMD) (Ragsdale et al., 2007): ${\displaystyle \sum _{n=1}^{\infty }=\left(\frac{xi-1+xi}{2}\right)\times \left(ti-ti-1\right)}$ where n is the total number of days over which sampling took place, x_i is the number of mealybugs counted on day i and t_i is the number of days since the initiation of sampling on day i.

Mealybug population build-up under each treatment was computed, and average CMD measures were compared between the respective treatments using a mixed modeling approach with plot as the random factor and time as the repeated measure. A mixed modeling approach was used to contrast fresh root yield (FRY), using treatment and variety as fixed factors. Plant survival rates were compared between treatments, using a Chi-square analysis. Where necessary, data were transformed to meet assumptions of normality and homoscedasticity, and all statistical analyses were conducted using SPSS. In cases where data could not be transformed to normality, non-parametric tests (e.g., Kruskal–Wallis) were used.

Country-wide yield changes

Crop production statistics were obtained through the Office of Agricultural Economics, Ministry of Agriculture & Cooperatives (Bangkok, Thailand). Yield measures were computed for 2006–2016, for a total of 51 cassava-growing provinces within Thailand, and annual weighted means were compared between successive years. Province-level yields were assigned different weights, based upon the local extent of cassava cultivation (in terms of harvested area) during a given year. Cassava crop yield can be impacted by agro-climatic conditions (e.g., temperature-related variables) and by attack of pests such as P. manihoti. To assess the impact of sustained A. lopezi releases from the 2011 cropping season onward, mean values of yields across all the cassava-growing provinces were regressed with explanatory variables which included rainfall, minimum and maximum temperature (obtained from Thai Meteorological Department, Bangkok, Thailand). Additionally, time (year for which yield observations and agro-climatic data were obtained) was also added as an explanatory variable to control for variation in the variables across years. In addition, a categorical variable representing the introduction of A. lopezi (dummy-coded as 1 for ‘presence’ for the 2011 and 2012 growing seasons, and dummy-coded as 0 for ‘absence’ for growing seasons 2008, 2009 and 2010) was also added as an additional explanatory variable in the regression model. Before proceeding for further analysis, the distribution of the response variable (i.e., yield) was tested and was identified to be normal (Shapiro test, p < 0.05). A step-wise regression approach (forward and backward) using a linear modeling approach was used to identify the model that best explains variation in yield. The model with the lowest Akaike information criterion (AIC) was selected. In the next step, the model with the lowest AIC score was compared with models containing interaction terms between time and other explanatory variables (i.e., temperature minimum, rainfall and A. lopezi introduction) both separately and simultaneously (see Table S1). Model with the least AIC score, and the highest adjusted R² value was selected and is described in the Results section. Regression analysis was performed in R (v 3.4.1) statistical computing environment. Additionally, R package “gvlma” was used to assess if the assumptions of regression were met by the selected model. Additional diagnostics of the selected model, such as determination of variance inflation factor (VIF) for detection of multicollinearity, the Non-constant Variance Score Test (i.e., test for heteroscedasticity of residuals over fitted values) was performed using R package “MASS” and “car”, respectively.

Multi-country pest & natural enemy survey

During continental-scale insect surveys from 2014 until 2017 (i.e., 5–8 years following the initial A. lopezi introduction), the mealybug complex on cassava largely comprised four non-native species: (1) P. manihoti; (2) the papaya mealybug Paracoccus marginatus Williams & Granara de Willink; (3) Pseudococcus jackbeardsleyi Gimpel & Miller; and (4) the striped mealybug Ferrisia virgata Cockerell. Phenacoccus manihoti was the most abundant and widespread mealybug species, and was reported from 37.0% (n = 549) and 100% fields (n = 52) in mainland Southeast Asia and Indonesia, respectively. Among sites, P. manihoti reached field-level incidence of 7.6 ± 15.9% (mean ± SD; i.e., proportion mealybug-affected tips) and abundance of 14.4 ± 31.0 insects per infested tip in mealybug-affected fields (or 5.2 ± 19.8 insects per tip across all fields) in mainland Southeast Asia, and incidence rates of 52.7 ± 30.9% and 42.5 ± 67.7 individuals per tip in Indonesia. Field-level incidence and population abundance were highly variable among settings and countries, reaching respective maxima of 100%, and 412.0 individuals per tip (Fig. 1).

When examining P. manihoti parasitism rates from a select set of sites, A. lopezi was present in 96.9% of mealybug-affected fields (n = 97) in mainland Southeast Asia, yet were only found in 27.5% sites (n = 40) across Indonesia. Among sites, highly variable parasitism rates were evident with dry-season rates of 16.3 ± 3.4% in coastal Vietnam, versus 52.9 ± 4.3% in intensified systems of Tay Ninh (also in Vietnam). In Indonesia, A. lopezi was found in 22.0% fields in Lombok (n = 9) and was absent from prime growing areas in Nusa Tenggara Timur (NTT). In sites where A. lopezi had successfully established, dry-season parasitism ranged from 0% to 97.4%, averaging 30.0 ± 24.0% (n = 110) (Fig. S1). In fields where A. lopezi had effectively established, mealybug pest pressure exhibited a negative regression with parasitism rate (F_1,98 = 13.162, p < 0.001; R² = 0.118).

Multi-year mealybug and parasitoid population assessment in Vietnam

Over the course of three years, we monitored P. manihoti abundance, field-level incidence and associated A. lopezi parasitism rates in Tay Ninh, southern Vietnam. Field-level incidence of P. manihoti ranged from 0% to 82%, averaging 24.8 ± 17.7% (mean ± SD) plants infested over two consecutive crop cycles. Mealybug incidence was significantly higher on older crops F_7,57 = 9.9; p < 0.0001 ), and rapidly increased during the dry season. Similarly, mealybug abundance (average 5.6 ± 5.0 individuals per tip) was higher during the dry season (F_1,63 = 9.10; P = 0.0037), and in crops older than six months, when compared to younger crops (F_7,57 = 2,694.06; P < 0.0001). From mid-2014 onward, mealybug populations remained at field-level abundance below 10–15 individuals per infested tip (Fig. 2). Furthermore, A. lopezi attained mean parasitism rates of 42.3 ±  21.7%, with maxima of 76.7 ± 28.9% during the early rainy season (Fig. 2). Overall, parasitism gradually increased over the dry season, up until crops were 4–6 months old.

Mealybug growth rates were significantly and negatively correlated with parasitism levels across the 8 sites studied (GLM w/ Normal error distribution and corrected for field: ${\chi }_1^2=125.4$; P = 0.0017; the field term was not significant; ${\chi }_7^2=0.2$; P = 1) (Fig. 3). The x-intercept of each per-field regression represents the parasitism level above which mealybug growth rates are negative and this value ranged between 0.38 and 0.69 for the 8 sites (average = 0.47  ± 0.09) (Fig. 3). Whilst A. lopezi was the sole primary parasitoid at this location, three hyperparasitoid species attacked it locally it at 2.79 ± 5.38% levels (as % of parasitized hosts).

Exclusion cage assays

Over the entire assay, P. manihoti populations under ‘full cage’ attained 48,318 ± 51,425 (n = 4; mean ± SD) and 7,256 ± 8,581 cumulative mealybug days (CMD) in ‘sham cage’ for one popular variety (i.e., R72) (Thancharoen et al., 2018; Fig. 4). For a second variety, KU50, P. manihoti attained 28,125 ± 32,456 CMD in a ‘full cage’ treatment, and 1,782 ± 1,073 CMD in ‘sham cage’. This compared to CMD measures in a ‘no cage’ control of 1,378  ± 1,039 and 342 ± 252, for R72 and KU50 respectively. CMD measures were significantly affected by treatment (F_3,189 = 240.752, p < 0.001) and time (F_6,189 = 113.347, p < 0.001), and the interaction term time x treatment (F_18,189 = 2.012, p = 0.011). Also, total CMD measures at the end of the trial significantly differed between treatments for both R72 and KU 50 (F_3,12 = 6.767, p = 0.006; F_3,12 = 11.152, p = 0.001, respectively).

Cassava yield parameters varied considerably under the four experimental treatments, and for both crop varieties (see Thancharoen et al., 2018). For Rayong 72, plant survival attained 37.5% under ‘full cage’ as compared to 75% and 87.5% under ‘no cage’ or ‘sham cage’ conditions, respectively (Chi square, χ² = 10.473, p = 0.015). Fresh root yield (FRY) was significantly affected by treatment (F_3,27 = 4.104, p = 0.016) and variety (F_1,27 = 4.364, p = 0.046). For R72 and KU50, FRY under ‘full cage’ was 74.6% or 71.2% lower than under ‘sham cage’ (Kruskal–Wallis, χ² = 8.344, p = 0.039; χ² = 19.134, p < 0.001, respectively), and respective yield reductions for both varieties were 77.2% and 67.8% compared to ‘no cage’ treatments.

Country-wide yield changes

During the 2009 dry season, P. manihoti attained its peak population in Thailand, with field-level incidence near 100% and abundance rates of hundreds of P. manihoti per plant on at least 230,000 ha (Rojanaridpiched et al., 2013). Over the subsequent 2009-10 cropping season, province-level crop yields dropped by 12.59 ± 9.78% nationwide (weighted mean: −18.2%) (Fig. 5). Furthermore, country-wide aggregate yields declined from 22.67 t/ha to 18.57 t/ha, and total production dropped by 26.86% to 22,005,740 tonnes of fresh root. Following the lowered crop output, prices for Thai cassava starch increased 2.38-fold at domestic prices in Thailand, and 2.62-fold at export prices (US$ FOB, Free On Board) (Fig. S2). However, in the following growing seasons (i.e., 2011 and 2012), there was a marked improvement in yield across all regions, as compared to previous years (Fig. 5). It was also from mid-2010 onward that country-wide mass releases of A. lopezi were carried out. To test if the A. lopezi introduction was responsible for this improvement in yield, and to differentiate P. manihoti-induced yield drops from agro-climatic climatic impacts and changes across years, regression analyses were carried out with yield as response variable and agro-climatic parameters, year (“time”), and A . lopezi presence as explanatory variables. Multiple regression analysis revealed that a model with interaction terms between time and all explanatory variables, i.e., time of introduction of A. lopezi and rainfall had the lowest AIC score and the higest explanatory power (adjusted R²) (Table S1). The model showed a significantly positive effect (F_7,183 = 8.641) of the interaction term Time ×Presence (i.e., ‘presence’ of A. lopezi and time , p < 0.01) on observed yields. The best model indicated that the introduction of A. lopezi (and not agro-climatic variables or changes in those variables over time) significantly increased yields across all cassava-growing regions during the 2011 and 2012 cropping seasons.