Mealybug species from Chilean agricultural landscapes and main factors influencing the genetic structure of Pseudococcus viburni

The present study aimed to characterize the distribution of mealybug species along Chilean agro-ecosystems and to determine the relative impact of host plant, management strategy, geography and micro-environment on shaping the distribution and genetic structure of the obscure mealybug Pseudococcus viburni. An extensive survey was completed using DNA barcoding methods to identify Chilean mealybugs to the species level. Moreover, a fine-scale study of Ps. viburni genetic diversity and population structure was carried out, genotyping 529 Ps. viburni individuals with 21 microsatellite markers. Samples from 16 localities were analyzed using Bayesian and spatially-explicit methods and the genetic dataset was confronted to host-plant, management and environmental data. Chilean crops were found to be infested by Ps. viburni, Pseudococcus meridionalis, Pseudococcus longispinus and Planococcus citri, with Ps. viburni and Ps. meridionalis showing contrasting distribution and host-plant preference patterns. Ps. viburni samples presented low genetic diversity levels but high genetic differentiation. While no significant genetic variance could be assigned to host-plant or management strategy, climate and geography were found to correlate significantly with genetic differentiation levels. The genetic characterization of Ps. viburni within Chile will contribute to future studies tracing back the origin and improving the management of this worldwide invader.

Scale insects are worldwide-distributed agricultural pests that cause major economic losses of billions of dollars every year 1,2 , either through direct impact on crops or quarantine export restrictions 3,4 . Commonly known as mealybugs, scale insects belonging to the family Pseudococcidae are particularly difficult to manage 5,6 . Their small size and cryptic habits render them very inconspicuous so that they easily escape phytosanitary inspections. For example, the obscure mealybug Pseudococcus viburni (Signoret) has expanded across the globe and is now present in over 60 countries 7 . Chemical control of pseudococcids is often ineffective due to their concealed nature and patchy distributions 4,8 . As a consequence, many management programs rely on classical biological control and especially the use of encyrtid parasitoids 6,9 . The response to particular biological control agents varies depending on the mealybug species 10 , so characterizing the diversity and distribution of mealybugs becomes important for both fundamental research as well as in applied pest management.
The spatial distribution of agro-ecosystems is a key factor determining the species composition and population genetic structure of agricultural pests 11,12 . Chilean highly structured agricultural landscape represents an excellent case study, with 83% of fruit orchards being found in a 700 km long stretch of land (between Coquimbo and Maule regions) and 20 fruit species concentrating about 90% of the cultivated area 13 . Geographical and environmental features of the country also contribute to further structuring Chilean agro-ecosystems. In the Norte Chico area (26° S to 32° S), which includes the administrative regions of Atacama (partly overlapping the Atacama Desert) and Coquimbo, there are semi-desert conditions with extreme climate variations. In the Central Chile area (32° S to 37° S), comprising the administrative regions of Valparaíso, Metropolitana, O'Higgins, Maule and the Northern part of Biobío, climate is mostly Mediterranean, with reliable alternation of warm, dry summers and cooler, rainier winter seasons 14 . The Norte Chico climate is mostly favorable to vineyards, while the Central Chile area includes both vineyards and fruit orchards (e.g. orange trees, apple trees, plum trees and pears). By running parallel to the Ocean, the coastal mountains and Andes cordillera further enclose vineyards and fruit orchards in isolated valleys, even within administrative regions.
National production of grapes and fruit orchards in Chile has significantly increased during the last decade 15 . Mealybug detection after phytosanitary inspections is responsible for over 45% fruit rejections in exports, and pest infestations must be controlled for the sustainability of the Chilean agro-industry 16 . Correct identification is essential for pest management, however previous reports on mealybugs have evidenced considerable difficulties when using morphology-based methods [17][18][19] . In a preliminary study using a DNA barcoding approach, Correa et al. 20 found Pseudococcus viburni (Signoret) to be the most common mealybug species in Chilean vineyards, followed by Pseudococcus meridionalis Prado and Pseudococcus cribata González. In contrast with previous studies though, no Pseudococcus longispinus (Targioni Tozzetti) or Planococcus ficus (Signoret) were found, which could result from the limited area being covered (two regions) and the fact that only vineyards were sampled 20 .
The present study aimed at characterizing the diversity and distribution of mealybug species infesting Chilean agro-ecosystems, including both vineyards and deciduous fruit crops. In order to achieve this goal, we completed an extensive survey between Atacama and the Northern part of Biobío (covering most cultivated areas in Chile) and using DNA barcoding methods to identify mealybug samples to the species level. Moreover, a fine-scale study of Ps. viburni genetic diversity and population connectivity was carried out using microsatellite markers. A comprehensive analysis of Ps. viburni samples allowed us to determine the relative impact of host plant, management strategy, geography and micro-environment on shaping the distribution and genetic structure of the pest.

Methods
Sampling and molecular identification. A total of 38 sampling sites were surveyed along Chile in 2012 and 2013 (Table 1). Sample collection was performed by taking one mealybug specimen per plant (20 -49 mealybugs per site) to avoid collecting individuals from a single family. Specimens were stored in absolute ethanol at − 20 °C until DNA extraction. Given that mealybug species are difficult to distinguish based on their morphology, DNA barcoding identification was performed. Genomic DNA was extracted from at least five individuals per locality using the DNeasy extraction kit (QIAGEN, Hilden, Germany), and the conserved 28S ribosomal marker was PCR-amplified using the primers S3690 (5′ -GAGAGTTMAASAGTACGTGAAAC-3′ ) and A4394 (5′ -TCGGARGGAACCAGC-TACTA-3′ ) 21 . PCR reactions were performed in a total reaction volume of 25 μ L and following the protocol of Abd-Rabou et al. 22 . The PCR mix was composed of 12.50 μ L 2X Tmix (QIAGEN, Hilden, Germany), 0.125 μ L of each primer (initial concentration of 100 μ M) and 10.25 μ L of ultrapure water. PCR conditions were: initial denaturation at 95°C for 15 min; 35 cycles of denaturation at 95 °C for 30 s, hybridization at 58 °C for 90 s, elongation for 60 s; final extension at 72 °C for 10 min. PCR products were then analyzed on a QIAxcel Advanced System (QIAGEN) and sent for sequencing to Genoscreen (Lille, France) or Beckman Genomics (Takeley, United Kingdom). PCR products were sequenced on both strands and consensus sequences and alignments were created manually with Bioedit version 7.01 23 .

Microsatellite Diversity Analyses.
Out of the initial 38 sampling sites, a total of 16 samples were identified as Ps. viburni (see results section) and kept for further analysis using the 21 polymorphic microsatellite markers developed by Correa et al. 24 . Genotyping was carried out in a total of 529 individuals from those Ps. viburni populations. Individual PCRs were performed in a total volume of 10 μ l, with cycling conditions as in Correa et al. 24 . PCR products were separated by electrophoresis using an ABI 3700 sequencer (Applied Biosystems) and 500 LIZ GeneScan TM size standard. Allelic profiles were obtained for each individual using the Genemarker TM v1.75 software (SoftGenetics LLC).
For each sampling locality, mean number of alleles, observed (H O ) and expected (H E ) heterozygosity were computed using the software GeneClass v2 25 . Allelic richness was estimated standardizing sample size with the R package standArich (available online at http://www.ccmar.ualg.pt/maree/software. php?soft= sarich). Deviations from Hardy-Weinberg Equilibrium (HWE) were tested for each locus and sampling locality using the software Genepop v4.2 26 . Per locus p-values were corrected for multiple testing within each population using the false discovery rate (FDR) procedure 27 . The presence of null alleles was evaluated with the FreeNA software, and 10,000 bootstrap replicates were used to test for significance 28 31 were performed using the software Genepop v4.2 26 . Levels of genetic differentiation between sample sites were also estimated using the DST statistic, which includes a multiplicative partitioning of diversity, based on the effective number of alleles rather than on the expected heterozygosity 32 . Pairwise DST statistic estimates and their significance using 1,000 bootstrap replicates were obtained with the pair.pops.Dest function in the DEMETICS package in R (http:// www.r-project.org/). Sampling localities were grouped with the neighbor-joining algorithm and the pairwise genetic chord distance 33  Genetic clusters were then inferred by minimizing deviation from HWE and linkage disequilibrium (LD) with the Bayesian clustering approach implemented in STRUCTURE 34 . Parameters were estimated assuming an admixture model with correlated allele frequencies and no prior location information. We carried out 20 replicate runs for each value of the number of clusters (K), set between 1 and 16 (i.e. the number of sampling sites). Each run consisted of a burn-in period of 200,000 Markov chain Monte Carlo (MCMC) iterations, followed by 1,000,000 MCMC iterations. The highest level of genetic structure was inferred through the Δ K approach as implemented in STRUCTURE Harvester 35,36 . Finally, a discriminant analysis of principal components (DAPC) was used to identify clusters of genetically-related genotypes. DAPC (like the sPCA method presented below) differs from the Bayesian method in STRUCTURE by not assuming HWE or linkage equilibrium 37 .
Spatial genetic structure analyses. Isolation by distance (IBD) was tested by performing a linear regression of the genetic distances (FST or DST) and the geographical distances (Km) and obtaining the Pearson's correlation coefficient. Matrices of genetic distance and geographical distance were subjected to a Mantel test with 10,000 permutations in Arlequin v.3.5 38 . A spatially-explicit multivariate method (sPCA: spatial analysis of principal components) was used to further explore the spatial patterns of genetic variability between sites 39 . Spatial connectivity networks were built using the sampling location data and three common algorithms (Delaunay triangulation, relative neighbor and sphere of influence) 40 . In order to select for the network that better fits the genetic data, the corrected Akaike Information Criteria (AICc) was computed for each graph using the ortho.AIC function of the spacemakeR package 41 . The connectivity graph with the highest AICc was used to test for the presence of spatial autocorrelation against the null hypothesis that allele frequencies are distributed at random. Spatial correlation can be either positive or negative, resulting in principal components with positive or negative eigenvalues. Principal components with significant positive eigenvalues indicate that sites are genetically more similar to their neighbors than expected by chance (global structures), while those with significant negative eigenvalues highlight local differentiation among neighboring sites 39 . Spatial autocorrelation (Moran's I) was tested for each component using the permutation procedure implemented in the spdep package and 1,000 permutations 41 .

Non-genetic factors influencing population structure. Several analyses of molecular variance
(AMOVA) were performed in order to evaluate the importance of non-genetic categorical factors on the distribution of genetic diversity. In particular, sample sites were grouped according to (i) administrative region (AT, CQ, VL, MT, OH or ML), (ii) type of agricultural management (conventional or organic production), and (iii) host plant (grape, apple, plum or pear). The significance for the AMOVAs was obtained through non-parametric permutation procedures (with 20,000 permutations) as implemented in Arlequin v3.5 38 . To examine whether non-genetic environmental factors may influence genetic differentiation levels in Ps. viburni, temperature and precipitation variables were recovered to 1-km spatial resolution from WorldClim version 1.4 (http://www.worldclim.org) 42 . The climate variables included temperature and precipitation parameters, namely: i) annual mean temperature, ii) maximum temperature of the warmest month, iii) minimum temperature of the coldest month, iv) annual precipitation and v) precipitation seasonality (Coefficient of Variation). A Spearman's rank correlation (ρ ) test was then carried out between pairwise genetic distances (FST or DST) and the Euclidean distance between climate variables for each pair of sampling sites. Spearman's correlation measures association based on the rank differences between two vectors and indicates how well a monotonic function describes their relationship.

Molecular identification of mealybug species. The present survey of scale insects infesting
Chilean crops identified four mealybug species (Planococcus citri, Pseudococcus longispinus, Pseudococcus meridionalis and Pseudococcus viburni). Ps. viburni was the most widespread pest, being found on vineyards, apple trees, plum trees or pear trees in six out of seven administrative regions sampled: Atacama (AT), Coquimbo (CQ), Valparaíso (VL), Metropolitana (MT), O'Higgins (OH) and Maule (ML) ( Table 1; Fig. 1). Ps. longispinus was also present in several regions and dominated apple fields south of Maule and in Biobío. Finally, Pl. citri and Ps. meridionalis were only found sporadically on vineyards or orange trees (Fig. 1).  (Table 2). After correction by the FDR procedure, seven sites still showed significant HWE deviations ( Table 2).  Table 3). All pairs of Ps. viburni samples showed significant differentiation levels after Fisher's exact tests (P < 0.05). The results obtained using the distance-based neighbor-joining approach showed a clear support for the separation of AT1 and CQ1 from the remaining samples (Fig. 2a). Despite most bootstrap values were low, a strong support was found for some groups sites, such as VL1-VL2, ML2-ML3, MT1-MT2 and VL3-VL4. Structure Harvester selected K = 3 as the best representative number of clusters based on likelihood and Δ K (Ln = − 16457; Δ K = 443.35). The first cluster included the two Northernmost sites (AT1 and CQ1), a second cluster was composed of samples from the Aconcagua (VL1 and VL2), Cachapoal (OH1), Colchagua (OH2 and OH3) and Curicó (ML1) valleys; and a third cluster was formed by the Casablanca (VL3 and VL4), Maipo (MT1 and MT2), Colchagua (OH5) and Curicó (ML2 and ML3) valleys (Figs 1,  2b). Finally, the DAPC method also support the presence of a cluster including AT1 and CQ1 samples, but indicated large overlapping between the samples assigned to the second (blue) and third (green) clusters given by STRUCTURE (Fig. 2c).

Microsatellite Diversity
Spatial genetic structure analyses. AMOVA results showed that a significant 7.73% of the total genetic variability can be explained by administrative regions (FCT = 0.074; P < 0.0001). A significant positive correlation was found between geographical distance and either FST (r 2 = 0.79; P < 0.0001) or DST values (r 2 = 0.75; P < 0.0001) when considering all the 16 sample sites, indicating the presence of IBD. When the IBD test was performed without the two most distant localities (AT1 and CQ1), a weaker but still significant pattern was also observed both using FST (r 2 = 0.068; P = 0.028) or DST values (r 2 = 0.121; P = 0.004).
According to the corrected AIC criterion, the best connectivity graph for explaining the Ps. viburni allele frequency data was obtained from the Sphere of Influence (SoI) algorithm (AICc = − 42.356; non-zero links = 34), which showed a lower AICc value than both the Delaunay triangulation (AICc = − 41.325;    Table 4). As for the effect of climatic environmental factors, Spearman's rank test indicated a strong correlation between genetic distances and differences in maximum temperature of the warmest month (ρ = 0.615; p-value < 0.001), minimum temperature of the coldest month (ρ = 0.645; p-value < 0.001), annual precipitation (ρ = 0.523; p-value < 0.001) and precipitation seasonality (ρ = 0.537; p-value < 0.001). However, the correlation between genetic distances and annual mean temperature was much smaller and non-significant after FDR correction (ρ = 0.239; p-value = 0.008).

Discussion
The distribution of insect pest species across agricultural ecosystems can be shaped by multiple factors, ranging from specific associations with the cultivated plant, environmental gradients or the use of pesticides. and was more abundant in Apple orchards south of that region. These different spatial distribution and host-plant preference patterns are in agreement with previous studies on the biology and systematics of Pseudococcidae 4 . Indeed, Ps. viburni belongs to the grape mealybug (Ps. maritimus) species complex, while the long-tailed mealybug (Ps. longispinus) and related species are commonly found on fruit orchards (e.g. avocado 43,44 and pear 45 ).
The host preference observed for different species did not extend to different populations within Ps. viburni, since no significant genetic variance could be assigned to host plant or management strategy. This lack of genetic structuring by host indicates that Ps. viburni populations are able to infest fruit orchards in an opportunistic way despite their preference for vines 4 . Rather than host plant or management strategy, our results highlight the importance of geography and environmental factors in shaping the intra-species population genetic structure of the obscure mealybug. The integrative analyses of micro-climatic variables and genetic markers show that the severity of climate (i.e. differences in the extremes of temperature or precipitation seasonality) is significantly correlated with genetic differentiation levels. Meteorological changes or edaphic factors have been previously proposed to modify the physiology of the plant and alter its resistance to scale insect attack 46,47 , but this is the first study to provide evidence on the importance of environmental factors on shaping scale insects distribution and genetic structure.
Climate differences might correlate with geographical distances, so that spurious genetic-climate correlations could result from an isolation by distance pattern. AMOVA, IBD and sPCA analyses all support the presence of reduced connectivity between Ps. viburni samples related to their spatial distribution. The fact that IBD remained significant even after excluding the two most distant populations indicates that genetic differentiation is also important at a local scale, where climatic differences are minimal. Therefore, geography cannot be discarded as a key factor shaping the population genetic structure of the obscure mealybug. In a recent study, the maritime pine bast scale (Matsucoccus feytaudi) was found to be highly structured geographically, and the authors relate this fact to the limited dispersal capacity of the insect and the patchy distribution of the obligate host 48 . Mealybug species are known to present low dispersal abilities, so long distance movements could only be driven by human activities and agricultural practices 4 . The peculiar shape of Chile, being only 177 km wide in average but 4300 km long, together with the low mobility of mealybugs, could well explain the observed geography-driven genetic structure of P. viburni.
The identification of well-differentiated groups within Ps. viburni is an important contribution to design effective control strategies for pest management both in Chile and elsewhere. Even though the exact center of origin for the obscure mealybug remains unclear, South America has been suggested as the most likely native area. Therefore, the genetic characterization of Ps. viburni within Chile will contribute to future studies tracing back the origin and improving the management of this worldwide invader (e.g. allowing to search for specific parasitoids within the source populations).