Unraveling the Genetic Structure of the Coconut Scale Insect Pest (Aspidiotus rigidus Reyne) Outbreak Populations in the Philippines.

Despite the fact that massive outbreaks of the coconut scale insect pest, Aspidiotus rigidus Reyne (Hemiptera: Diaspididae) are inflicting significant economic losses to the Philippines' coconut industry, little is known about the population and dispersal history of this invasive pest in the country. Here, we examined the genetic diversity, structure and demographic history of A. rigidus sampled from localities with reported outbreaks from 2014 to 2017. We analyzed the genetic structure of the outbreak populations using mitochondrial COI and nuclear EF-1α markers. Both markers and all methods of population genetic structure analyses indicate clear differentiation among the A. rigidus populations separating the north from the southern regions of the Philippines. Very low or no genetic differentiation was observed within and amongst the populations per geographic region indicating two unrelated outbreak events of the pest originating from two genetically uniform populations isolated in each respective region. Historical data supports the resurgence of an established A. rigidus population in the south which could have been driven by sudden climatic changes or human-induced habitat imbalance. Our study provides valuable information on the genetic differentiation of the two A. rigidus groups that would be useful for developing and implementing biological control strategies against this pest in the Philippines.

1. Introduction

Insect pest outbreaks are characterized by an explosive increase in an insect population occurring over a relatively short period [1]. Large and rapid alterations in the environment or changes in the intrinsic genetic or physiological properties of individual organisms within a population can result to the resurgence of insect pests to outbreak-level status [2,3]. Likewise, insect outbreaks may occur when non-native species have no natural enemies or a few inefficient natural enemies, and if the local beneficial species are unable to suppress them in the area of introduction [4,5]. The invasive success of pest species may be determined by both the biology and environmental factors promoting its spread in a suitable area [6,7]. A better understanding of the source population, route and the mechanism of spread could provide valuable insights for designing and implementing quarantine strategies to understand the invasion success and decline of outbreak populations [4,8].

In 2009, the Philippines suffered a devastating coconut scale insect outbreak damaging the coconut palms on the provinces of Luzon (northern region of the Philippines) and currently on some areas in Mindanao (southern region) inflicting significant economic losses to the country’s coconut industry. The diaspidid scale insect Aspidiotus rigidus Reyne (Hemiptera: Diaspididae) resides on the underside of the leaf, blocks the stomata and sucks plant sap strongly reducing the plant’s photosynthetic activity leading to a characteristic yellowing and drying of the leaves. Severely infested coconut palms dry up and die within six months or less [9]. Prior to the renewed interest on A. rigidus due to the outbreak in the Philippines, historical and observational data on the spread of the invasive coconut pest has been scarce. Other than the most recent observation published by Watson et al. [10] the last known study on the biology of this invasive species was conducted by Reyne [9,11] with full documentation of the outbreak on the island of Sangi (North Celebes) in Indonesia from mid-1925 to 1928. The outbreak recorded by Reyne [9] came to an end naturally after two years due to reduced female fecundity and high mortality of the immature stages. The decrease in the A. rigidus population may have been associated with natural enemies that, over time, regulated the pest population. However, the outbreaks from the localities infested with A. rigidus in the Philippines took longer to recover [10], e.g., six years for the northern province of Batangas, and are still on-going for the southern areas, i.e., Basilan, Zamboanga Peninsula, and the Caraga Region.

The introduction of A. rigidus to the Philippines, and its spread was believed to be either by wind or accidental transportation of infested plants, coconut planting materials and products [10]. Infestation in the northern provinces of the Philippines spread rapidly from its initial local report in Tanauan, Batangas in the Calabarzon Region (Luzon) in 2009, reaching nearby coconut planted areas throughout the region. These outbreaks lasted for at least three years [10] and were reported to be manageable by 2015 [12]. The more recent outbreak in the southern region, specifically in Basilan, started in early 2013 [10], implying its direct connection with the northern outbreak. However, given the means of spread via wind wherein crawlers are dispersed from one area to another [10], it is highly improbable that the infestation from the northern region reached the infested southern islands meanwhile moving past other provinces planted with coconut palms en route. In addition, transport of infested plants from the northern provinces was highly unlikely given the national attention focused on quarantine and management strategies against the spread of the coconut scale insect during the outbreak [12,13].

It is also likely that A. rigidus has been in the country as a minor pest and regulated by natural enemies. Based on historical reports, Lever [14] reported sightings of A. rigidus in the Philippines as long ago as 1969, and Velasquez [15] recounted that the pest was probably highly confined in the southern part of the country. It was more likely that the source of the observed A. rigidus came from the island of Sangi in Indonesia given its relative closeness to Mindanao. In time, the immigrant A. rigidus could have established a resident population complemented by natural enemies limiting its colonization outside the area of introduction. Changes in anthropogenic, biotic interactions or climatic factors can influence a population’s rise to an outbreak level [3,16]. The recent outbreak observed in the southern part of the Philippines may have been caused by a sudden rise in the abundance of the supposedly established A. rigidus population due to factors such as human-induced habitat imbalance e.g., excessive use of pesticide affecting the natural enemies controlling the pest population, or climate change such as a prolonged dry spell, which may induce changes in the local biotic community.

Inference of the source population, route and the mechanism of spread of A. rigidus in the Philippines needs further assessment and confirmation. Tracing the history of an invasion or identifying the geographic origin of a pest population can be done by characterizing population-level genetic variation using molecular markers e.g., [17,18,19,20]. Sequencing selected gene fragments, e.g., mitochondrial COI is a traditional population genetic tool providing insights on dispersal pathways and population structure. The mitochondrial cytochrome oxidase (mtCOI) gene and the nuclear protein-encoding gene elongation factor 1α (EF-1α) have been commonly used in studies investigating the origin [21,22], or inference of phylogenetic relationships [23,24] of various invasive diaspidid species.

Here, we aimed to assess the population genetic structure and demography of the outbreak populations of the coconut scale insect, A. rigidus in the Philippines. Given the historical documentation of the pest in the southern region and the relatively extensive and rapid spread, but faster recovery, of the infestations in the northern region compared to the southern outbreaks i.e., Basilan and Zamboanga Peninsula, we hypothesize the presence of two distinct genetic groups for the outbreak events isolated within each geographic region. A population genetics approach is a useful tool to examine whether the northern and southern coconut scale insect outbreaks originated from immigrant or resident populations. To test this hypothesis, we utilized sequences of the mtCOI gene and the nuclear protein-encoding gene - EF-1α to investigate the genetic structure and diversity of A. rigidus populations from localities with documented outbreak-level infestations in the Philippines from 2014 to 2017. Furthermore, we employed a coalescent genealogy approach to provide additional evidence on the demographic relationship of the outbreak A. rigidus populations between the northern and southern geographic regions in the Philippines.

2. Materials and Methods

2.1. Sample Collection

Aspidiotus rigidus populations were sampled at seven localities with reported coconut scale insect outbreak across the Philippines from 2014 to 2017 (Figure 1). The northern localities sampled were Orani, Bataan (BT; N14.769786, E120.454510), Nagcarlan (NG; N14.158930, E121.413670) and San Pablo (SP; N14.056420, E121.333300), Laguna, Tanauan (TN; N14.098870, E121.091330) and Talisay (TL; N14.093340, E121.010730) Batangas. The southern localities are Basilan (BS; N06.707853, E121.983358) and Zamboanga (ZB; N06.993166, E121.927963). Detailed information regarding location and sample collection is collected in Table 1.

Figure 1

Map of the seven localities with reported A. rigidus outbreaks in the Philippines from 2014 to 2017. The insect rearing facility of the Biological Control Research Unit of De La Salle University labeled “AR”. Dots indicate sampling locations. Northern localities: Orani, Bataan (BT), Nagcarlan (NG) and San Pablo (SP), Laguna, Tanauan (TN) and Talisay (TL), Batangas; Southern localities: Basilan (BS) and Zamboanga (ZB). See Table 1 for the more detailed information regarding location and sample collection information.

Table 1

Sampling localities of the outbreak A. rigidus populations. N, number of individuals with mtCOI and EF-1α sequences; H, haplotypes indicated in Figure 3.

Locality		Code	Collection Date	mtCOI		EF-1α
				N	H	N	H
Northern Region	DLSU-STC, Laguna a	AR	July 2017	35	H1	8	H1
	Orani, Bataan	BT	September 2015	21	H1	0	—
	Nagcarlan, Laguna	NG	January 2015	24	H1	11	H1
	San Pablo, Laguna	SP	December 2014	21	H1	8	H1
	Talisay, Batangas	TL	December 2014	35	H1	18	H1
	Tanauan, Batangas	TN	December 2014	29	H1	8	H1
Southern Region	Isabela, Basilan	BS	November 2016	44	H2	16	H2; H3
	Zamboanga City, Zamboanga	ZB	April 2017	96	H2	6	H4

aA. rigidus reared on mangosteen at the DLSU-STC BCRU rearing facility from samples collected on the outbreak population in Orani, Bataan.

Mature female scale insects were identified as A. rigidus based on the characteristic distribution of egg skins, which for this species occurs along the posterior or pygidial half of the insect body (Figure 2). Non-parasitized adult females were carefully selected from infested leaves and preserved in 95% ethanol before molecular analysis. To further confirm identification as well as the purity of samples, A. rigidus collected from Orani, Bataan were reared at the rearing facility of the Biological Control Research Unit (BCRU) located at De La Salle University (DLSU), Science and Technology Complex (Binan City, Laguna) on mangosteen (Garcinia mangostana L.), a differential host of A. rigidus observed not to support multiple generations of A. destructor. A phylogenetic analysis was employed (see the molecular analysis section below) to confirm the identifications of the field-collected samples by comparing them to A. destructor and mangosteen-reared A. rigidus sequences.

Figure 2

Representative adult female A. rigidus from different outbreak areas: (A) Southern Tagalog Region (Laguna, Cavite, and Batangas); (B) Orani, Bataan; (C) Basilan; and (D) Zamboanga City. The arrows point to the egg skins, which for this species is characteristically distributed along the posterior or pygidial half of the insect body.

2.2. DNA Extraction, PCR Amplification, and Sequencing

Genomic DNA was extracted individually using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) following the manufacturer’s guidelines. Extraction was performed by crushing the insect body of each individual in separate microcentrifuge tubes using a micropestle. DNA concentration and quality were assessed by spectrophotometry (NanoDrop 2000 spectrophotometer, ThermoScientific, Wilmington, DE, USA). The mitochondrial COI gene was amplified using the forward primer PcoF1 designed for scale insects by Park et al. [25] and the standard reverse primer LepR1. The nuclear gene EF-1α was amplified using the forward primer EF-1α of Morse and Normark [26] paired with the EF2 reverse primer [27]. PCR reactions were performed in a 25 μL reaction mixture containing 10× Buffer, 2.5 mM dNTP mixture, 25 mM MgCl2, 10 pmol of each primer, 1 U of Taq DNA polymerase (TaKaRa Bio Inc., Japan), and 2 to 50 ng of template DNA. PCR thermocycling was performed in a T100™ Thermal Cycler (Bio-Rad, CA, USA). Following the conditions from Park et al. [28], the mtCOI gene was amplified with an initial denaturation step at 95 °C for 5 min, followed by 5 cycles of 94 °C for 40 s, annealing at 45 °C for 40 s, extension at 72 °C for 1 min and 10 s, and another 35 cycles of denaturation at 94 °C for 40 s, annealing at 51 °C for 40 s, extension at 72 °C for 1 min and 10 s, and a 5 min final extension at 72 °C after the last cycle. Furthermore, after an initial denaturation at 95 °C, with a denaturation at 95 °C for 30 s and extension at 72 °C for 2 min every cycle, a touch-down procedure was performed for the amplification of the EF-1α gene following the protocol of Morse and Normark [26] in which the initial annealing temperature of 58 °C was decreased by 2 °C every three cycles until a final temperature of 42 °C was reached, then held for 18 cycles followed by a 5 min final extension at 72 °C. PCR products were visualized in 1.5% agarose gels stained with Midori Green Direct (Nippon Genetics Co. Ltd., city, Japan), and cleaned using the QIAquick PCR Purification Kit (Qiagen, Tokyo, Japan). Samples were sent to Eurofins Genomics (Eurofins Genomics Co., Ltd., city, Japan) for Sanger sequencing to produce both forward and reverse fragments.

2.3. Genetic Diversity and Population Structure

Sequences were assembled using CodonCode Aligner v. 5.1.5 (CodonCode Corporation, Dedham, MA, USA). Concatenated sequences were deposited in GenBank with the following accession numbers: MK862684-MK862988 for mtCOI and MN603810-MN603884 for EF-1α before the subsequent molecular analysis, the sequences were aligned via MAFFT v. 7.409 [29], and the ambiguously aligned regions were excluded using GBlocks 0.91b [30]. Sequence polymorphisms for both mtCOI and EF-1α gene were assessed. The number of variable sites (S) and haplotypes (h), average number of nucleotide difference (k), haplotype diversity (Hd), and nucleotide diversity (Pi) of the two marker genes were calculated in DnaSP v. 6.10.04 [31]. The hierarchical analysis of molecular variance (AMOVA) was implemented in Arlequin v. 3.5.2.2 [32]. Two geographic groups were defined: The northern (Luzon) and the southern (Mindanao) groups prior to the analysis. Population pairwise F was computed in Arlequin v. 3.5.2.2 using 1000 permutations.

2.4. Demographic Inference

Tajima’s D and Fu’s Fs statistic tests were estimated to infer demographic history and dynamics in each population, for the two geographic groups, and for all populations grouped together in Arlequin v. 3.5.2.2 for both mtCOI and EF-1α datasets. Also, Fu and Li’s D* and F* test statistics were computed in DnaSP v. 5.0 to determine departures from the mutation-drift equilibrium [33]. Parameters of demographic expansion included moment estimators of time to the expansion Tau, effective population size before expansion (Theta0, θ), and effective population size after expansion (Theta1, θ) between the observed and expected mismatches. The adjustment to a model of population expansion was estimated from the sum of squared deviation (SSD) and the raggedness index (r) in Arlequin v. 3.5.2.2.

2.5. Gene Flow Analysis and Median-joining Networks of Haplotypes

To test migration history between the two geographic groups, we calculated Bayes factors from the marginal likelihoods estimated in MIGRATE v. 4.4.0 [34,35] based on both mtCOI and EF-1α datasets. MIGRATE-N utilizes marginal likelihoods to compare and order structured population models [35]. The program provides estimates of historic gene-flow with the assumption that populations have reached mutation-migration-drift equilibrium. We tested eight possible models of migration history. Model 1 allows migration between the two groups, with the populations assumed to exist for a very long time. Model 2 presents a migration from northern to southern group, while model 3 presents the opposite situation and with the populations assumed to exist for a very long time. Model 4 assumed one panmictic population encompassing the northern and southern groups. Model 5 allows divergence among populations within southern group splitting from the northern group, and migration from north to south, with the northern group existing for a long time and the southern group recently splitting off. Model 6 is a mirror image of model 5. Model 7 is similar to model 5 except that no interaction occurred between the two groups after the split. Model 8 represented the opposite situation of model 7. Similar parameters were used to run all models.

A Bayesian search strategy was performed with the following parameters: one long chain (10,000 trees) with a burn-in of 5000 iterations. A static heating scheme with 4 chains was applied using temperature parameters set by default with a swapping interval of one. Bayes factors were calculated via “BF” implemented in carlopacioni/mtraceR, a package for analyzing migrate-n outputs in R v. 3.5.2. Log Bayes factors (LBF) of all models were calculated by comparing against the model that has the highest log-likelihood. The models are ranked based on LBF and calculated model probability.

Median-joining (MJ) networks [36] for the two markers were constructed to estimate the genealogical relationship in A. rigidus haplotypes via PopART v. 1.7 [37].

2.6. Phylogenetic Analysis

Since no site variation was observed between the sequences in each population (except for the EF-1α sequences from Basilan), six representative samples per population, a total of 48 mtCOI sequences along with 42 plus all 16 Basilan EF-1α sequences were chosen for the phylogenetic analysis. Sequences of A. destructor collected from coconut palms, identified based on the circular distribution of egg skins were selected as an outgroup (accession number: MK862989 for mtCOI and MN603885 for EF-1α). The Akaike information criterion corrected for sample size (AICc) was implemented to find the best fitting evolutionary model for phylogenetic reconstruction via jModelTest v. 2.1.10 [38]. The evolutionary model for the mtCOI sequences was TIM2+G, while TrNef was the model for the EF-1α sequences. Maximum likelihood (ML) tree inference was performed in RAxML-NG v. 0.5.1 [39] with 1000 bootstrap replicates.

3. Results

3.1. Genetic Diversity and Population Structure

All samples identified based on the characteristic distribution of egg skins were confirmed as A. rigidus by comparing the sequences with the mangosteen-reared A. rigidus, and A. destructor sequences. DNA sequence analysis of all the concatenated 647-bp mtCOI sequences of 305 individuals from seven A. rigidus outbreak populations collected from 2014 to 2017 in the Philippines, with the mangosteen-reared samples revealed only two distinct haplotypes (h), separated by 31 polymorphic sites (s) (Figure 3). Haplotype diversity (Hd) was calculated to be 0.50 +/− 0.005 SD. Average number of nucleotide difference (k) was 15.5 and nucleotide diversity (Pi) was 0.024 +/− 0.00025 SD. No variation was found in the sequences of samples collected per populations. Reared A. rigidus and samples collected from the five populations of northern group, i.e., Orani, Bataan (BT), Nagcarlan (NG) and San Pablo (SP), Laguna, Tanauan (TN) and Talisay (TL), Batangas are grouped into one haplotype. Samples from the two populations of southern group, i.e., Basilan (BS) and Zamboanga (ZB) were grouped together in the second mtCOI haplotype (Table 2). For the nuclear EF-1α gene, 75 concatenated sequences of 1007-bp length generated 14 polymorphic sites (s), with four haplotypes (h). Like the mtCOI sequences, all samples from the northern group clustered into one haplotype. For the southern group, samples from ZB grouped into one haplotype, while BS were separated into two haplotypes. Haplotype diversity (Hd) was calculated to be 0.48 +/− 0.064 SD. Average number of nucleotide difference (k) was 5.7 and nucleotide diversity (Pi) was 0.0057 +/− 0.00061 SD. Except for BS, all other localities have no sequence variation per locality.

Figure 3

Median-joining network of the A. rigidus populations from 305 individuals for the mtCOI gene (A), and 75 individuals for the protein-coding EF-1α gene (B), showing location and frequency of haplotypes. Each circle represents an observed haplotype; circle size indicates the number of individuals observed; the colors correspond to sampling localities. The total number of mutations, Eta presented as hatch marks.

Table 2

Parameters of genetic diversity and demographic analysis of the two population groups.

Gene	Group	N	S	h	Haplotype a	k	Hd (SD)	Pi (SD)	D b	Fs	D * b	F *
mtCOI	Northern	165	—	—	H1	—	—	—	—	—	—	—
	Southern	140	—	—	H2	—	—	—	—	—	—	—
	All	305	31	2	—	15.45	0.4980 (0.005)	0.0239 (0.00025)	5.85 ***	59.59	2.04 **	4.43
EF-1α	Northern	53	—	—	H1	—	—	—	—	—	—	—
	Southern	22	2	3	H2; H3; H4	0.87	0.6880 (0.039)	0.0009 (0.00010)	1.33	0.99	0.85	1.13
	All	75	14	4	—	5.71	0.4770 (0.064)	0.0057 (0.00061)	2.83 **	12.81	1.55 *	2.38

a Haplotype data by DnaSP v. 6.10.04; b Parameters with statistical test: * indicates p < 0.05; ** indicates p < 0.02; *** indicates p < 0.01.

In addition, genetic diversity parameters were calculated per geographic group. The mtCOI sequences for both groups, and the EF-1α sequence of the northern group showed no sequence variation. The EF-1α sequences from the southern group had three haplotypes with an estimated Hd of 0.688 +/− 0.039 SD, with a k value of 0.087 and a low Pi value of 0.00086 +/− 0.00010 (Table 2). Both the median-joining haplotype network and ML inferred trees present distinct two and four haplotypes for the mtCOI and EF-1α dataset, respectively (Figure 3 and Figure S1).

AMOVA analysis indicated a highly structured genetic variability of 100% and 97.63% variations among the groups for mtCOI and EF-1α dataset, respectively. There were zero, or a relatively small percentage of variation among populations within groups and within populations. Except for the source of variation among populations within groups in the mtCOI data, AMOVA showed that significant genetic structure occurred in A. rigidus at various hierarchical levels (Table 3). Pairwise F values varied from zero to 1.00 for mtCOI, and zero to 0.98 for the EF-1α dataset. The differentiation between populations was only significant when the comparison was between a northern and a southern population (Table S1). Moreover, pairwise F values between the two groups revealed a high value of 1.00 and 0.97 for the mtCOI and EF-1α dataset, respectively. As shown in the ML trees (Figure S1), the phylogenetic analyses of both markers were consistent with the results of the analyses above. Samples clustered according to their geographic group, with the BS and ZB EF-1α sequences in three separate nodes.

Table 3

Partitioning of genetic variation at different hierarchical levels.

Gene	Source of Variation	df.	Sum ofSquares	VarianceComponents	Percentageof Variation	Fixation Indices
mtCOI	Among groups	1	2347.869	15.50000 Va	100	FCT = 1.00000 *
	Among populations within groups	6	0	0.00000 Vb	0	FSC = 0.00000
	Within populations	297	0	0.00000 Vc	0	FST = 1.00000 **
EF-1α	Among groups	1	202.010	6.45405 Va	97.63	FCT = 0.97630 *
	Among populations within groups	5	5.199	0.09876 Vb	1.49	FSC = 0.63040 **
	Within populations	68	3.938	0.05790 Vc	0.88	FST = 0.99124 **

* indicates p < 0.05; ** indicates p < 0.01.

3.2. Demographic History and Gene Flow

For both datasets, neutrality test computation for all samples showed positive values for Tajima’s D, Fu and Li’s D*, Fu’s Fs, and Fu and Li’s F*, and were significant for the first two parameters (Table 2). Estimations per population for these parameters were mostly zero or positive but not significant, suggesting neither population expansion nor purifying selection in these populations. Estimations of the SSD and r parameters both returned zero values, except for the EF-1α sequences from Basilan, with a significant SSD of 0.0283 (p < 0.001) and a non-significant r of 0.2871. Other demographic parameters such as Tau, θ and θ index, are presented in Table S2.

Results of the analysis in migrate-n were presented in Table S3. We found contrasting results for the two markers employed using Bayes factors to compare the eight models of dispersal. For the mtCOI dataset, model 7 was ranked best with a probability of 0.996. For the EF-1α dataset, model 3 was ranked best with a probability of 1.000.

4. Discussion

The main aim of the present study was to elucidate the genetic structure and demography of the coconut scale insect pest A. rigidus from selected localities in the Philippines with reported heavy infestations collected from 2014 to 2017. Both the mtCOI and EF-1α markers and all methods of population structure analyses revealed strong differentiation among the A. rigidus populations separating the northern (Luzon) outbreak from the southern (Mindanao) region. The separation of the populations by geographic region and the observed lack of genetic variability within populations were represented graphically in the median-joining network and phylogenetic analysis employed in the study (Table 1 and Figure 3).

4.1. Genetic Structure of A. rigidus: Evidence of Coconut Scale Insect “Superclones” in the Philippines

Our results indicate the existence of two mitochondrial, and four nuclear haplotypes (one northern and three southern). Genetic population clusters result from multiple source populations contributing to an insect pest outbreak [8]. However, we only observed two clusters separating the outbreak populations into their respective geographic regions. Also, genetic variation was either absent or very low within and amongst the populations of the northern and the southern region, implying that populations from each region consisted of a single genotype. Hence, the presence of two distinct A. rigidus single genotype populations or “superclones” [40] in the Philippines which supports our hypothesis on the occurrence of two genetically unrelated outbreak events in the country.

Armoured diaspid scale insects were observed to have unusual genetic systems such as diverse types of paternal genome elimination and parthenogenesis [22]. Specifically, several aspidiotine insects have obligate [24,41] or cryptic parthenogenetic populations as observed in Aspidiotus nerii Bouche [21,22]. Accordingly, A. rigidus was observed to reproduce parthenogenetically [9,10]. Yellow winged adult males are seen in outbreak populations but the sex ratio varies widely with males thought to be non-functional [9]. Parthenogenetic reproduction has been thought to be the leading driver to the dominance of “superclones” across space and time in other insects such as aphids (Hemiptera: Aphididae) [40,42]. As exemplified with our findings, some invasive insect pests have been found to depend on clonal population structures to successfully invade and multiply in a broad range of niches e.g., [43]. A highly specialized clonal genotype of a strictly asexual population of the pea aphid, Acyrthosiphon pisum Harris was the main reason influencing the demographic success of the pest in central Chile [44]. Cifuentes et al. [45] found no genetic variation within populations of the tomato leaf miner, Tuta absoluta Meyrick (Lepidoptera: Gelechiidae) spreading through South America and ultimately reaching the Mediterranean Basin, with only a single genetic type. Likewise, a well-established invasive population of the oleander aphid, Aphis nerii Fonscolombe were reported as having extremely low genetic diversity in the southern United States, with a “superclone” population supposedly obligatorily asexual [46]. Caron et al. [47] reported two widespread, invasive and strictly parthenogenetic “superclones” of the sawfly, Nematus oligospilus Forster (Hymenoptera; Tenthredinidae) dominating willows in three countries in the southern hemisphere i.e., South Africa, New Zealand, and Australia.

4.2. Demographic History: Resurgence of Resident Population or Recent Introduction?

Based on historical reports, Lever [14] claimed that the more invasive coconut scale insect A. destructor rigidus (now A. rigidus) reported by Reyne [9] in Indonesia was also present in the Philippines. However, Velasquez [15] did not observe its occurrence across the archipelago and reported a highly probable confinement of the pest in the southern region of the country. Given this historical evidence, we assume that the southern populations have existed for a long time. This scenario might have supported our dispersal model for the mitochondrial sequences, except that the first reported sighting of A. rigidus in Tanauan, Batangas, Luzon was in 2009 [10] with no historical evidence of resident populations in the past. This suggests that the northern populations were most probably a recent introduction event from a different source.

A Bayesian search strategy was performed to assess the migration history between the northern and southern populations. However, our results from the mitochondrial and nuclear datasets are difficult to reconcile. We reiterate that in this historic gene-flow analysis, the two groups were assumed to have reached mutation-migration-drift equilibrium. Despite the contrasting results, both models indicate that sequences from the two groups do not belong to one panmictic population. Even so, given the difference in the divergence or migration pattern of the models for each marker, we were unable to unequivocally infer the source of each outbreak population. Methodological assumptions [48] in the program Migrate-n, just like other coalescent-based approaches, did not take into account another source of migrants, or that ancestral variation may come from populations that were not considered in the analysis. Hence, the inference of the possible source of the northern outbreak population needs further exploration.

On the other hand, lower genetic variation is expected for younger populations due to founder effects and genetic bottlenecks during colonization and establishment [49]. Invasive or recently introduced species have been reported to exhibit reduced genetic variation e.g., [50,51]. Introduced populations are usually small so decreased genetic diversity is expected, and are often less variable than the source population which contributes to the invasive success of the species [45]. The nuclear marker revealed the existence of three southern haplotypes, with samples from Basilan having two distinct haplotypes. Genetic variation amongst the populations was very low and the Zamboanga EF-1α sequences are differentiated compared with that from Basilan with few nucleotide substitutions (Figure 3B).

Genetic variation was already low amongst the southern populations for the nuclear DNA, but in comparison to the uniformly genetic northern population (Table 2 and Table 3), our findings infer that the southern A. rigidus population was relatively older than the northern one. Alongside previous historical reports, the level of genetic variation between the geographic regions supports our hypothesis of an existing resident A. rigidus population in the southern part of the Philippines. Local insect populations have the potential to outbreak due to anthropogenic and environmental changes [1,3]. Similar observations on insect pests have been reported in the literature. A notable example by Kobayashi et al. [8] reveals that the multiple nationwide outbreaks of the native populations of the mirid bug, Stenotus rubrovittatus Matsumura (Heteroptera: Miridae) in Japan were induced by changes in the agro-ecosystem without invasion of populations from other areas. Populations of the pest were also genetically isolated by distance and separated into genetic clusters occupying spatially segregated regions. Additionally, temporal fluctuations of pest insects in agroecosystems could be driven by various factors [2]. Pesticide application may induce the resurgence of native pest insect populations by reducing the abundance of natural enemies or by the removal of competitive species in the area e.g., [52,53]. Weather conditions can also trigger insect outbreaks due to the dramatic changes in pest abundance. Ward and Aukema [54] reported that the cyclic outbreaks of the native tree-killing bark beetle, Dendroctonus simplex LeConte (Coleoptera: Curculionidae) on Tamarack (Larix laricina (Du Roi) K. Koch) in Minnesota, USA are climate-driven, being specifically associated with warmer and dryer years, more likely in areas with prior defoliation. Schwartzberg et al. [55] simulated climate warming and observed warming-induced phenological shifts in the forest tent caterpillar, Malacosoma disstria Hübner (Lepidoptera: Lasiocampidae) in relation to the phenology of its host trees. These findings illustrate the mechanisms by which anthropogenic and climatic changes induce outbreaks from native insect pests.

4.3. Future Prospects: Applications for Biological Control

The information regarding the amount of genetic variation in pest species could be useful because insects with reduced genetic variation have been reported to be more easily controlled than those with more genetic diversity [56]. For instance, populations with multiple resistant phenotypes are more difficult to control via pathogenic biological control agents compared to populations lacking diversity [57]. In 2014, Almarinez et al. [58] discovered a new encyrtid species, Comperiella calauanica Barrion, Almarinez, Amalin (Hymenoptera: Encyrtidae) parasitizing A. rigidus infesting the coconut palms in a forest-based resort in Calauan, Laguna (Luzon, Philippines) with no history of chemical insecticide applications. The parasitoid has been observed to be host-specific, with up to 80% parasitisation rate. The use of C. calauanica for biological control implementations through conservation and augmentation in the field have been conducted in the Philippines [58,59]. Utilizing the natural variation of the pest insect, or lack thereof, would be vital in the development and implementation of biological control strategies for integrated pest management programmes against this coconut scale insect pest.

5. Conclusions

The current opinion about the origin of the coconut scale insect outbreak in the Philippines was a recent introduction of A. rigidus from other countries of native range and spread via wind dispersal or importation of infested planting material from the northern region to the south given the timeline of the outbreak reports. However, our results strongly suggest the separation of two distinct groups, the northern and southern A. rigidus from the outbreak populations collected from 2014 to 2017 in the Philippines. Very low or no genetic differentiation was observed within and amongst the populations per geographic region, indicating two unrelated outbreak events of the pest species originating from two genetically uniform or “superclone” populations currently isolated in each respective region. Historical data supports our assumption on the current resurgence of an established A. rigidus population in the south. Given the lack of historical information supporting the existence of an established A. rigidus populations in the northern region, we disregard the possible resurgence of a persisting native population and suggest that the outbreak possibly resulted from a recent introduction of a non-native population. Assessment of the possible source population of the northern outbreaks needs further exploration.

The use of mtCOI and the nuclear EF-1α markers showed no or very low genetic differentiation for all A. rigidus populations. Other robust and more informative genetic markers such as microsatellites could provide further genetic information during study of the invasive coconut scale population. Further studies should also include more expansive sampling, taking into consideration other possible sources of A. rigidus such as Indonesia [10] and Vietnam [24]. This would provide a more robust and stringent population and gene flow estimation of A. rigidus in the Philippines. Nevertheless, our findings provided an initial important genetic basis and information for designing and implementing biological control strategies against the invasive coconut scale insect pest A. rigidus in the Philippines.