Tracing the dispersal route of the invasive Japanese beetle Popillia japonica

doi:10.21203/rs.3.rs-1497043/v1

Download PDF

Research Article

Tracing the dispersal route of the invasive Japanese beetle Popillia japonica

https://doi.org/10.21203/rs.3.rs-1497043/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

The Japanese beetle, Popillia japonica, is a highly polyphagous Scarabeidae native to Japan that colonized North America and Azores in the last century and has recently invaded Italy and Switzerland. Considering its economic impact to the horticulture and turfgrass industries, this species was ranked within the EU priority pests list in 2019.

In order to reconstruct the source of introductions of this pest, we investigated the genetic variability of P. japonica in its native and invaded areas worldwide by analyzing 9 microsatellite loci and two mitochondrial genes, COI and CytB. There are two populations of Japanese beetles, from a limited area within central Japan, that are likely the source of the North American population. Moreover, the microsatellite data and resulting phylogeographic reconstruction suggests that the two European populations originated from independent introductions from southeast (Azores) and northeastern (Italy and Switzerland) North America ₍₁₎.

(1) For simplicity, in this paper North America refers to Canada and USA

Different Popillia japonica populations were identified in native area.
The North American Popillia japonica population originated from a limited area on Honshu between Gunma, Ibaraki and Tochigi Prefectures.
The Azorean and European mainland populations of Japanese beetle could be considered independent introductions from northeast and southeastern North America respectively.

The Japanese beetle, Popillia japonica Newman, is an invasive insect belonging to the family Scarabaeidae. Native to Japan, it is considered one of the most widespread and destructive pests of turfgrass and ornamental plants worldwide (Potter and Held 2002).

Japanese beetle is a highly polyphagous species that can feed on more than 300 species of wild and cultivated plants (EFSA et al. 2019) as a larva and adult. Due to the economic impact to the horticulture and turfgrass industries, this species was ranked within the EU priority pests list in 2019 and it is now receiving a lot of attention from many European National Plant Protection Organizations (EU Commission Delegated Regulation 2019).

At mid-latitudes, P. japonica displays a univoltine life cycle, spending most of its life cycle as a larva in the soil (Potter and Held 2002). In the summer, females lay eggs in the soil, the larvae develop through 3 instars and emerge as adults at the beginning of the following summer (EFSA et al. 2019). The flight period typically peaks in mid-July and adult activity terminates by the end of October (EFSA et al. 2019). Female P. japonica are highly fecund, typically ovipositing about 40–60 eggs during her lifespan (Fleming 1972). Furthermore, P. japonica is not only a strong flier, able to disperse about 3–24 km per year (Fox 1932), but it can also relocate passively to uninfested areas through hitchhiking on non-host commodities or vehicles (EPPO 2016).

In Japan, P. japonica occurs at high population densities in northern Honshu and Hokkaido. It is also common in Kyushu, Shikoku and southern Honshu (Fleming 1972; Klein 2008). Locally, it was not considered a pest until golf courses were damaged because of larval activity. The first report of P. japonica outside of Japan occurred in 1916 when it was first detected in southern New Jersey, USA (Fleming 1976). In the following years, it colonized the entire eastern USA and progressively expanded to the Midwest. With the exception of Florida, which probably has an unsuitable climate, the Japanese beetle became established in all states east of

the Mississippi River by 1998 (Potter and Held 2002). Only sporadic occurrences are currently reported from western USA, where several eradication campaigns have been carried out using chemical insecticides (Clair and Kramer 1989; Hungate et al. 2016). In Canada P. japonica was first detected in 1938 and is now present in the eastern provinces of Ontario, Quebec, Nova Scotia, New Brunswick, and Prince Edward Island (CFIA 2009). In 2017, an infestation was discovered in Vancouver, British Columbia but this population is currently being eradicated (EPPO 2018). Outside of North America, P. japonica reached the Azores (Portugal) in the '70s and is now established on several islands in this region, although at different population densities (Potter and Held 2002). The first outbreak in mainland Europe occurred in Ticino Natural Regional Park, northern Italy, in 2014 (Pavesi 2014). This area is located between the Piedmont and Lombardy regions in an area where two airports, one commercial and one military, are present. In 2017, P. japonica was also reported from southern Switzerland, possibly as an expansion of the Italian infestation (EPPO 2017). Scattered reports of P. japonica in the Korean peninsula, China and India may on the other hand have arisen from misidentifications (Potter and Held 2002).

Losses due to P. japonica feeding on wild and cultivated plants can be substantial. Although it is uncertain how much economic damage Japanese beetle causes in Europe, it is estimated that more than 450 million dollars are allocated annually to manage the spread of this pest and to replace damaged plants in the US (Potter and Held 2002). In Italy, several measures are being employed to limit the distribution of Japanese beetle. Natural and commercial entomopathogenic nematodes and fungi have been tested and applied to infested soils (Benvenuti et al. 2019; Marianelli et al. 2018; Mazza et al. 2017; Paoli et al. 2017; Torrini et al. 2020) to control larvae. Attract and kill devices, made of long-lasting insecticide-treated nets with a pheromone lure, along with “attract and infect” devices equipped with entomopathogenic fungi and also a pheromone lure, have been developed and deployed in infested areas to control adult populations (Marianelli et al. 2019; Regione Piemonte 2021). Even with these control tactics, P. japonica has invaded as much as 15,000 square kilometers in northern Italy and in the neighboring Swiss Ticino Valley since 2014 (Mori et al. 2022).

The experience of the last decades suggests that knowledge of the genetic structure and the routes of introduction of an invasive species can assist in developing effective monitoring programs and limit additional introductions (Nardi et al. 2005; Roderick and Navajas 2003). Moreover, novel approaches are being developed that use molecular targets for the genetic control of pests, further underlining the importance of genetic data and data on the geographic structuring of genetic variability (Scott et al. 2020).

In order to determine how the current distribution of Japanese beetle occurred, our research investigated the genetic variability of P. japonica in its native and invaded areas worldwide by analyzing 9 microsatellite loci and two mitochondrial genes, COI and CytB.

Sampling and DNA extraction

Adult and larval P. japonica were collected from several areas in its current global distribution. More than 215 individuals were collected and the sampling coordinates, date of collection, and name of collector were recorded for each specimen (Table S1). Following the manufacturer’s recommended protocol, DNA was extracted from adult thoracic muscles and the last three larval metamers, using Blood and Tissue DNA Mini Kit (QIAGEN). The final elution step was performed in 100 µl of AE buffer supplied in the commercial kit.

Microsatellite development

Low coverage sequencing of the P. japonica genome was performed at BioFab Research s.r.l. (Rome, Italy) to produce a total of 2.31 million contigs (from 78bp to 978Kbp) that encompass 1.57Gbp of DNA sequences (data not shown). Automatic extraction of microsatellite loci was performed by this lab to identify 1,620 regions that include sequences (di-, tri- and tetra-nucleotide) longer than 10 repetitions. Based on the initial observation that different contigs share regions of sequence similarity and that ancillary repeats, simple sequences and homopolymer runs are frequently observed in flanking regions, candidate loci were subject to a filtering step before primer design.

Blastn ver. 2.6.0+, function makeblastdb; (Camacho et al. 2009) was used to compare the 1,620 regions (repeat and 200bp flanking sequences on each side) against the genomic assembly under default settings that, in terms of stringency, roughly correspond to a match of 100bp with 90% similarity. Regions that produced more than one match against the genomic assembly were discarded.

The occurrence of sequences characterized by low complexity in the flanking regions, that could complicate primer design and sizing, was evaluated using a modified Lempel-Ziv algorithm (Orlov and Potapov 2004) (http://wwwmgs.bionet.nsc.ru/mgs/programs/lzcomposer/index.html). This measure, developed for the study of non-coding genomic regions, is stable for complex sequences, but displays a significant drop in the presence of dispersed repeats, palindromes, homopolymer runs and other sequences of low complexity. Sequence complexity was evaluated in a window of 40bp with 10bp increments over 200bp flanking regions on each side. Regions with at least one segment of 40bp with associated complexity lower than 15 were discarded. The threshold was adjusted based on visual observation of the results. Furthermore, only sequences characterized by perfect di- and tri-nucleotide repeats, recorded in contigs of length > 3,000bp, with an A + T content < 70% and devoid of homopolymer runs > 5bp in the 200bp flanking the repeat were retained. The resulting 71 regions were visually inspected and 20 were selected for primer design and testing. Primers were designed using the primer-blast (Ye et al. 2012) implementation of Primer3 (Untergasser et al. 2012).

Optimization of candidate loci

Twenty candidate primer pairs were used for testing and optimization of amplification conditions. Amplifications were performed using a GoTaq G2 Flexi DNA Polymerase (Promega, Madison, WI, USA) following the manufacturer’s specifications, with cycle 95°C 2m; 95°C 1m, 50–58°C 1m, 10s 72°C 1m 30s x 35 cycles; 72°C 10m, 4°C 10m on a GeneAmp PCR System 2700 (Applied Biosystems, Foster City, CA, USA). Following optimization of PCR conditions, fragments that produced clear amplicons were cloned in a pGEM-T Easy Vector System (Promega, Madison, WI, USA) and sequenced to confirm locus specificity. Ten loci that produced optimal amplifications were selected for genotyping. Microsatellite sequences have been deposited in GenBank under accessions OL802935/OL802944

Genotyping

Genotyping was performed on 16 individuals for each of the 7 sampling areas (112 individuals in total). All forward primers (BioFab Research s.r.l., Rome, Italy) were marked with 6-fam dye. Amplifications were performed as above, with annealing temperature as indicated in Table 1. Genotyping was performed at BioFab Research s.r.l. (Rome, Italy) on ABI3730 DNA Analyzer (ThermoFisher Scientific, Waltham, MA, USA) using a 1:200 dilution of PCR products. The LIZ conjugated ladder Gene Scan 500 (ThermoFisher Scientific, Waltham, MAS, USA) was used for sizing, with the exclusion of the 250bp band following the manufacturer’s instructions. Scoring of peaks was performed using Geneious (v.11.1.3, Biomatters, Auckland, New Zealand). The region between ladder bands 50-500bp was evaluated and peaks were called using the 3rd order least squares approximation. Bins were automatically inferred. Automatically detected peaks were manually revised to account for locus specific idiosyncrasies (e.g. tendency to form stagger bands) and read into a spreadsheet. One locus (F) was excluded because the frequent presence of smearing instead of amplification bands. The final dataset is available as Supplementary Data.

Microsatellite analysis

Basic locus statistics were calculated using GenePop v. 4.7 (Raymond and Rousset 2001; Rousset 2008). Clustering of individuals into putative populations was obtained in Structure v. 2.3.4 (Falush et al. 2003; Pritchard et al. 2000) using the admixture and correlated model, as well as the locprior option (Hubisz et al. 2009) as suggested in the software documentation for small sized datasets.

Fifteen independent runs were conducted for each possible value of K (1 to number of sampling areas in different analyses) with a burn-in of 50.000 and 100.000 MCMC generations after burn-in. Optimal K was selected using the Evanno method (Evanno et al. 2005) as implemented in Structure Harvester v. 0.6.94 (Earl and vonHoldt 2012). Among multiple replicates for a given K, the one with the highest likelihood was used for visualization in Structure Plot v. 2.0 (Ramasamy et al. 2014). The analysis was initially conducted on the full dataset. In order to deploy knowledge on the sequences of introduction/establishments in different areas, subsets of the data pertinent to specific phases/events of introduction were also analyzed separately (Japan only, Japan + USA, USA + Azores + Italy).

Amplification of mitochondrial targets, sequencing and data analysis

Fragments of mitochondrial genes COX I and CytB were amplified using the LCO1490/HCO2198 primers (Folmer et al. 1994) and CP1/CB-N-11367 primers (Harry et al. 1998; Simon et al. 1994), respectively (Table 1). Amplifications were performed in a total volume of 50 µl containing 25 µl of HS DreamTaq Master mix (2X) (Thermo Scientific) 0.6 mM of each primer and 10 ng of template DNA. The thermal protocol used for the amplifications was 94°C 3m (for HS activation); 95°C 30s, 40–45°C 30s, 72°C 45s x 45 cycles; 72°C 5m. PCR products were checked on 2,0% agarose gel, purified with QIAquick PCR purification kit (QIAGEN), and sequenced in both directions at BioFab Research s.r.l. (Rome, Italy). Sequences were assembled using Geneious (v.11.1.3, Biomatters, Auckland, New Zealand). Samples with incomplete or low-quality sequences in one or both mitochondrial loci were discarded in the subsequent analyses.

Alignments of the resulting 187 sequences were performed with ClustalOmega and manually inspected. Haplotype diversity (Hd), the average number of nucleotide differences (k), nucleotide diversity (π), Tajima’s D (Tajima 1989), and Fu’s Fs (Fu 1997) neutrality indices were calculated using Arlequin v. 3.5.2.2 (Excoffier and Lischer 2010) for both loci.

Maximum likelihood (ML) analyses were performed on individual loci and the concatenated dataset to evaluate the congruence between the two loci. The appropriate substitution matrices were chosen based on JModelTest v. 2.1.10 considering the corrected Akaike information criterion. Trees were computed using MEGA-X software adopting TN93 as substitution matrix in all three cases (COI, CytB, and concatenated loci) and a moderate branch swapping algorithm. A bootstrap test with 1000 pseudoreplicates was conducted.

Haplotype networks were reconstructed using the TCS algorithm (Clement et al. 2002) in PopArt (Leigh and Bryant 2015) for three COX I–CytB combined datasets: all samples, Japanese samples and those from introduced areas.

The Bayesian reconstruction of phylogeographical relationships between samples was performed using *BEAST2 v. 2.6.4 software (Bouckaert et al. 2014) plus bModelTest v. 1.2.1 (Bouckaert and Drummond 2017), model-selection and GEO_SPHERE v. 1.3.0 (Bouckaert 2016) as additional packages. Data were arranged in three partitions, two containing nucleotide sequences and the third containing geographical coordinates of the sampling area. Samples were arranged in 39 different OTU named from the geographical regions where specimens come from. The most suitable substitution matrix for each nucleotide partition was TN93 and the coalescent exponential growth was applied as a demographic tree prior to a strict molecular clock model. MCMC was run for one billion generations, sampling every 5,000 generations and discarding 10% as burn-in.

Dataset evaluation and basic statistics

From the initial 215 samples, we obtained successful amplification and complete sequencing of both mitochondrial loci in 183 Popillia japonica, 3 Popillia quadriguttata and 1 Popillia mutans specimens. Sequences were submitted to GenBank; accession numbers are reported in Table S1.

The analyzed regions of COX I and CytB mitochondrial genes of Japanese beetle were 657 and 744 nucleotides long respectively and did not include any gaps or missing data. Genetic indices and neutrality tests calculated for both loci are summarized in Table 2. Tajima’s D and Fu’s Fs values on COX I and CytB loci were negative (only COX I Tajima’s D was not significant) confirming the recent Japanese beetle expansion.

Genotyping produced clearly identifiable allele data for 9 out of the 10 microsatellite loci initially targeted, however there was a fair amount of data missing (only 5 loci were obtained in Italian specimens). Peak identification was straightforward and required minimal manual curation from the automated Geneious scoring. Loci were generally quite variable, with 6 to 29 alleles per locus overall and values of expected heterozygosity (He) from 0.34 to 0.91 (Table 1). The Hardy Weinberg exact test, conducted by locus and by population, rejected equilibrium (p < 0.01) in all but three instances (i.e. loci A, B, I).

Microsatellite data

Analysis of the full dataset (7 sampling areas) suggested K = 2 as the optimal number of clusters, with a secondary signal at K = 4. If two clusters are assumed (Fig. 1a), southern Japan is occupied almost exclusively by one cluster (green in Fig. 1a) while east and central-west Japan appear to be characterized by a mixture of both clusters in similar proportions. All invading populations belong to the latter cluster (orange), with a minimal reminiscence of the former in southeast USA. Partitioning in four clusters led to the identification of a) a first group (green) almost exclusive to southern Japan; b) intermixing of three groups in east and central-west Japan (orange, azure, purple), with a dominance of the former in central-west Japan; c) intermixing of the latter two groups (azure, purple, not orange) in the USA and their separation in the Azores and Italy, that appear populated by one different group each.

Analysis of the Japanese dataset (3 sampling areas) suggested K = 2 as the optimal number of clusters. One cluster corresponds unequivocally to southern Japan and the other to east and central-west Japan (Supplementary Fig. S2b).

Analysis of the Japan + USA dataset (5 sampling areas) suggested K = 2 as the optimal number of clusters, with a secondary signal at K = 4. Assuming 2 clusters, only one is identified in southern Japan while a combination between this and a second cluster characterizes east and central-west Japan. Only the latter is present in the USA (Supplementary Fig. S2c). Assuming 4 clusters, most sampling areas appear to carry an intermix of groups, with no clear discernible pattern. Of the 4 groups present in Japan, only 2 were introduced to the USA, and these are absent from southern Japan (Supplementary Figure S2d).

Analysis of the USA + Azores + Italy dataset (4 sampling areas) suggested K = 2 as the optimal number of clusters, with a secondary signal at K = 4. Assuming 2 clusters, these appear to be amalgamated in the USA but are separated in Azores and Italy, where a different cluster is observed in each location (Supplementary Fig. S2e). Assuming 4 clusters, each of the 4 sampling areas is characterized by the dominance of one cluster, with a higher level of intermixing in the USA than in the Azores and Italy. While differences in the less frequent groups are limited and not readily discernible, the cluster characterizing the Azores is more frequent in the US southeast than in the northeast, while the one characterizing Italy is more frequent in northeastern USA (Supplementary Fig. S2f).

Mitochondrial data

The genetic trees based on the two mitochondrial genes COX I and CytB obtained using ML were congruent and did not have topological conflicts, but displayed weak support for the regional subdivision of Japanese samples and therefore, do not have sufficient resolution to determine the topology of samples from non-native areas (Supplementary Figures S3 and S4). The TCS network constructed with COX I – CytB combined data obtained from the entire distribution of P. japonica included 84 nodes, but only 65 were vouchered haplotypes, with an average distance between nodes of 2.34 units with moda 1, and only a few nodes located very far from the main core of the network. The majority of haplotypes (n = 53) occur in Japan, but only 3 of these are also present in the colonized areas while 10 haplotypes only occur outside the native area (Fig. 2; see Tab S1 for full location information).

Two separate TCS networks using haplotypes found in samples collected in native and non-native areas respectively were computed to refer the haplotypic data to geographical regions of origin. The TCS network obtained only with Japanese samples connects 53 vouchered haplotypes with an average distance of 2.79 units (Fig. 3). Haplotype groups are geographically restricted to specific Japanese regions, e. g.: group IX was found only in the southern Japanese islands of Kyushu and Tsushima, while groups I and III in the center of Honshu. Other groups have a wider distribution such as IV and VII that were observed from the center of Honshu to Hokkaido. Although some groups were widespread across Japan, no single haplotype or group was observed throughout the archipelago.

The TCS network obtained from the introduced area displayed less variability compared to the Japanese one (n = 13), with only 8 haplotypes observed exclusively in North America (USA and Canada), two haplotypes among Azorean and Italian beetles, and one in the Swiss population (Fig. 4). Only three North American haplotypes were found in Japan (Hap_34, Hap_38, and Hap_48) where they were located in a geographic region that comprises the Gunma, Ibaraki, and Tochigi prefectures. The most common and widespread haplotype in the North American continent was Hap_38 while Hap_48 was found only in a limited region in USA between New Hampshire and Maine. Only 2 haplotypes were identified in Azorean Japanese beetles analyzed (Hap_38 and Hap_59). The same haplotypes of Azores were identified in the Italian outbreak of P. japonica, although with substantially different frequencies, while only Hap_59 was found in samples from Switzerland.

The phylogeographical reconstruction based on multipartition data (COX I, CytB, and geographical coordinates of sampling sites) obtained by Bayesian inference supports the subdivisions previously found between mitochondrial haplotypes and geographical coordinates. Moreover, the Bayesian tree was able to partially resolve the topology of non-native areas (Fig. 5a). Haplotypic groups previously identified in the native-area network were generally monophyletic, except for Group V that was definitively paraphyletic, and Groups II and III that were not well-resolved at a posterior probability threshold of 0.85. Specimens collected outside the native range occurred in two clades that are substantially superimposable to Group II and Group III. In North America, specimens belonging to Group III were found in the US Midwest, eastern US and eastern Canada, while Group II specimens were collected from the US Midwest, eastern US and western North America. The recent introductions in British Columbia, Oregon, Idaho and Colorado belong to Group II.

Popillia japonica has been a fairly neglected phytophagous species until it came to the limelight following the invasion of North America and the Azores in the last century, and Italy more recently. A rich scientific literature has since accumulated on outbreaks and their spread, as well as on the phenology of the species in different ecological settings. And a consensus of the human-mediated introductions from its native range in Japan to the North America and from here, to the Azores and mainland Europe has also occurred. Nevertheless, a reconstruction of the historical invasion routes of this species using a phylogeographic approach had been missing until now.

Genetic analyses of Japanese samples of Popillia japonica identified 2 different lineages in the native area, the first in Kyushu and Tsushima and the second across Shikoku, Honshu and Hokkaido. This separation between south and north-central Japan was jointly supported by microsatellites and by mitochondrial data, with the identification of two different genetic clusters in the two areas. Haplotype group IX, characterizing the southern islands, is separated by 13 substitutions from Group VIII (Shikoku) and 33 substitutions from the core of the network (Honshu and Hokkaido). This differentiation may have been promoted by the presence of the Seto Inland Sea, which could have acted as a natural barrier to gene flow separating Kyushu from Shikoku and Honshu over historical times. A parallel differentiation was, on the other hand, not observed between Shikoku and Honshu, with samples clustering in the same group based on microsatellite data and extensive sharing of haplotypes belonging to Group V across the two islands. Whereas the entire Honshu was recovered as a single assemblage based on microsatellite data, mitochondrial haplotypes recovered additional structure in the area. The highest haplotypic richness was observed in central Honshu, suggesting that this may be the origin of this species, with other Japanese islands having been subsequently colonized through a process of range expansion that has occurred at some point in the past. Haplotypic data from samples collected in the northern part of Japan (Hokkaido and Akita Prefecture in Honshu) were in turn characterized by Group IV and Group VII haplotypes, further supporting and old expansion process from the central region of Honshu to Hokkaido. Taken together, these observations suggest an old process of differentiation of P. japonica in its native area, with an origin in Honshu and secondary differentiation in northern and southern islands.

In the light of the biological features and the possibility for dispersal of P. japonica, its presence in North America is undoubtedly associated with human activities and trades that took place in the last century. Microsatellite analysis identified the origin of the North American outbreak from the east and central-west Japan populations. Based on mitochondrial data, the presumptive geographic origin of the invading population could be narrowed down further to central-eastern Honshu. Of the 13 haplotypes found in North America, only 3 were identified in Japan, namely Hap_34, Hap_38, and Hap_48, with the first two being the most frequent overall in North America and having spread and differentiated to produce a series of ancillary haplotypes. Other haplotypes occur sporadically in different sampling sites, with frequencies under 2%. It is therefore possible to hypothesize that Hap_34 and Hap_38 characterized the original introduction, with other satellite haplotypes being generated by casual mutation during the colonization process. These two most frequent North American haplotypes were identified only in a limited geographical region of Japan, in the center of Honshu island between the Ibaraki, Tochigi, and Gunma Prefectures. As such, this geographical region can be considered as the most likely origin of the invasive specimens associated with the first North American outbreak of P. japonica.

Being present in Japan in the same region (Nagano Prefecture), Hap_48 may be part of the original introduction, having followed the same route hypothesized for Hap_34 and Hap_38. Nevertheless, the observation that it is found at low frequency in North America, has not differentiated locally to produce satellite haplotypes and its geographical distribution is restricted to Maine (USA), gives credence to the alternative possibility that it may be the outcome of an independent outbreak that occurred more recently in a limited geographical area in the US.

Microsatellite analysis of the Azorean and Italian samples allowed only a weak definition of the presumptive area where European outbreaks originated. Adopting a 4-cluster model, it was possible to suggest that Azorean Popillia may have originated from southeastern USA and the Italian population from northeastern USA. Independent confirmation of the former hypothesis was obtained by analyzing mitochondrial data. Of the two haplotypes found in the Azorean Japanese beetles, Hap_38 is common in the US, while Hap_59 is observed in Colorado, New Jersey, and Canada's British Columbia. This latter differs from the more common Hap_38 by only 2 substitutions, with the intermediate haplotype Hap_62 identified in samples from North Carolina (USA). A possible origin of P. Japonica outbreak in the Azores could therefore be tracked to the US east coast. In fact, following from the idea of a continuum spread during the North America colonization, the first invaders characterized by Hap34 and Hap_38 could have spread in the south where Hap_62 and the subsequent Hap_59 were generated. It is therefore possible that some insects were accidentally introduced to the Azores more than 50 years ago from the coastal southeast of the US. This model would also explain the low haplotypic variability observed in the Azores as a consequence of a double founder effect.

The Italian outbreak is characterized by the same haplotypes of the Azorean population, nevertheless their haplotypic frequencies are totally different. Furthermore, the Italian microsatellite pattern seems to be more related to the northeast region of North America rather than the southeast. As such, and even allowing some shift in genetic frequencies associated with low population number during the invasion, it is unlikely that the Italian outbreak originated from the Azores. Of the three locations where Hap_59 was identified in North America, and in line with the indication of microsatellite analysis, New Jersey appears the more likely source, as Colorado and British Columbia are recent detections (2018 and 2017 respectively) and the pest is not considered established locally (CFIA 2020; EPPO 2018).

Circumstantial evidences suggest a possible correlation between planes and P. japonica, in particular the position of the initial Italian outbreak in an area between Milano Malpensa International Airport and Cameri military airport. Moreover, Hap_59 and Hap_38 haplotypes are present in sympatry only in New Jersey, where the Newark International Airport occurs, historically the largest airport in the US and currently part of the largest airport system word-wide. Altogether these observations allow us to hypothesize a possible origin of the Italian specimens from the New Jersey area mediated through accidental transport by plane.

The Swiss outbreak is monomorphic for Hap_59, which is the dominant haplotype in Italy. This, together with the geographical continuity between the two sites, strongly suggests that the spread of Italian P. japonica in the neighboring Ticino area of Switzerland was the origin of the Swiss outbreak.

In order to add the information of distances between sampling areas, we performed a Bayesian phylogeographical reconstruction using both sampling geographical coordinates and COX I and CytB mitochondrial data. Microsatellite data were not included in this analysis due to the long and intensive computation required (an average of 20 days for 1 million generations). In the native range the previously identified groups were monophyletic with the exception of Group V whose paraphyly can be due to a spotted geographical distribution (Fig. 5b). The outbreaks discovered in the western region of North America were originated only by specimens belonging to Group II, as a result of introductions from the east to west. While the eastern Canadian populations only belong to Group III and arose by establishment of populations from southern locations. Even if it is possible a different migration patterns occurred between the two north American populations of P. japonica, this evidence could also be also explained by the existence of two other human-mediated spreading pathways of Japanese beetle through USA and Canada, but still unidentified.

The genetic data reported here, including rejection of HW equilibrium in the majority of microsatellites loci and negative values of Tajima’D and Fu’S estimators that indicate an excess of rare alleles in mitochondrial haplotypes, suggest, in genetic terms, a situation of non-equilibrium. This is in line with multiple evidence, including recent phytosanitary reports (EPPO 2021), suggesting that Japanese beetle is in a demographic expansion phase, leading to new outbreaks and spread to neighbouring areas. However other factors, besides demographic expansion, could be at play to determine allelic disequilibrium, such as the strong founder effect that can take place in outbreaks and the population bottleneck caused by human activities aimed at reducing the impact of this pest. Altogether, this suggests some caution in the interpretation of genetic dissimilarities based on frequency data, as these are subject to change quickly and randomly during such processes. The shared presence of specific haplotypes that display a restricted distribution is, on the other hand a more reliable signature of shared ancestry, as it is less affected by genetic disequilibrium.

The Japanese beetle is an important pest, enlisted in the alert list in several countries around the world due to its polyphagy, possibility to create major damage to cultivated and wild plants and proclivity for long range successful introductions. In this work, we have described an old differentiation of the species in its native area with sufficient resolution to track back the origin of the North American P. japonica outbreak to a limited region in the center of Honshu. Our reconstruction found some support for a distinct origin of European P. japonica, with the Azorean outbreak presumably originating from southeastern USA and the Italian-Swiss outbreak possibly related to an independent introduction from northeastern USA. Our results also support the hypothesis that airports are high-risk entry points for new accidental introduction of this species. As such, the improvement of methods to evaluate the presence of Popillia adults and larvae in traded plants and commodities may be warranted.

Acknowledgements

Authors would thank all collectors for their precious contribution: Katsura Ito, Isao Takenaka, Shin-ya Ohba, Haruki Tatsuta, Ken Funayama, Yoshihiro Kunitomo, Kouichi Oyamada, Ryusuke Ishikawa, Kiichi Shimizu, Shingo Tanaka, Takuma Takanashi, Satoshi Kitabayashi, Koji Noge, Yoriko Saeki, Graham S. Thurston, Allison Barr, Ronald Neville, Jody, Green, Alan Eaton, Bill Jones, Claudio Gratton, Diana N. Kearns, David Held, Johnson Donald Steinkraus, Frank Krell, Sarah Zukoff, Paulo Castrovillo, Albrecht Koppenhöfer, R. C. McDonald, Charlene Donahue, Doug Terhune, Nadeer N. Youssef, Paulo Borges, Cristina Marazzi, Chun Lin Li and Chen Ri-Zhao. Authors would also thank Ilaria Cutino and Chiara Leo for their support.

Funding

This work has received funding from the regional project “COBIPO—Biological Control of Popillia japonica” [Piedmont Region (D n. 910 – 2015/12/29)] and European Union’s Horizon 2020 research and innovation programme “IPM POPILLIA” Integrated pest Management of the invasive Japanese Beetle, Popillia japonica, under grant agreement No. 861852.

Author Contributions

AS conceived experiments, AS, FP, FN, KS, TK, AC, LM planned sampling, AS, NF, II made laboratory work, AS and FN analyzed data, AS, FP, NF, II, TK, AC, LM wrote the manuscript and GB, PFR, LM provided found for the research, TK revised English language. All authors read and approved the manuscript.

Conflict of interest: Authors declare that they have no conflict of interest.

Ethical approval: This study does not contain any experiments using any animal species that requires ethical approval.

Consent to participate: Consent was given by all participants included in the study.

Consent for publication: All authors consent to the publication of this manuscript, figures, tables and supporting materials in Journal of Pest Science.

Data availability: All sequence data were available in GenBank database, microsatellite data were available in Supporting materials.

Benvenuti C et al (2019) A new device for auto-disseminating entomopathogenic fungi against Popillia japonica: a study case. Bull Insectol 72:219–225
Bouckaert R (2016) Phylogeography by diffusion on a sphere: whole world phylogeography. PeerJ 4:e2406. https://doi.org/10.7717/peerj.2406
Bouckaert R et al (2014) BEAST 2: A Software Platform for Bayesian Evolutionary Analysis. PLoS Comput Biol 10:e1003537. https://doi.org/10.1371/journal.pcbi.1003537
Bouckaert RR, Drummond AJ (2017) bModelTest: Bayesian phylogenetic site model averaging and model comparison. BMC Evol Biol 17:42. https://doi.org/10.1186/s12862-017-0890-6
Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL (2009) BLAST+: architecture and applications. BMC Bioinformatics 10:421. https://doi.org/10.1186/1471-2105-10-421
CFIA (2009) Questions and answers: Japanese beetle (Popillia japonica). http://www.inspection.gc.ca/english/plaveg/pestrava/popjap/queste.shtml. Accessed 4 November 2021
CFIA (2020) Appendix 1: Regulatory status of areas in Canada and the United States for Japanese beetle (Popillia japonica). http://www.inspection.gc.ca/plants/plant-pests-invasive-species/directives/horticulture/d-96-15/appendix-1/eng/1346826626609/1346826990603
Clair DJ, Kramer VL (1989) Japanese beetle. In: Dahlsten DL, Garcia R (eds) Eradication of exotic pests: analysis with case histories. Yale University Press, New Haven, CT, pp 89–107
Clement M, Snell Q, Walker P, Posada D, Crandall K TCS: estimating gene genealogies. In: Parallel and Distributed Processing Symposium, International, 2002. IEEE Computer Society, pp 0184–0184
Earl DA, vonHoldt BM (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour 4:359–361. https://doi.org/10.1007/s12686-011-9548-7
Schrader EFSA, Camilleri G, Ciubotaru M, Diakaki RM, Vos M S (2019) Pest survey card on Popillia japonica. EFSA Supporting Publications 16:1568E. https://doi.org/10.2903/sp.efsa.2019.EN-1568
EPPO (2016) PM 9/21(1) Popillia japonica: procedures for official control. EPPO Bull 46:543–555. https://doi.org/10.1111/epp.12345
EPPO (2017) First report of Popillia japonica in Switzerland. EPPO Reporting Service Num. article: 09-160
EPPO (2018) Popillia japonica found in Vancouver (British Columbia, Canada). EPPO Reporting Service Num. article: 05–092
EPPO (2021) EPPO datasheets on pests recommended for regulation: Popillia japonica. https://gd.eppo.int
EU Commission Delegated Regulation (2019) 1702 of 1 August 2019 supplementing Regulation (EU) 2016/2031 of the European Parliament and of the Council by establishing the list of priority pests. Official journal of the European Union L. 260/8, 11/10/2019
Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software structure: a simulation study. Mol Ecol 14:2611–2620. https://doi.org/10.1111/j.1365-294X.2005.02553.x
Excoffier L, Lischer HEL (2010) Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour 10:564–567. https://doi.org/10.1111/j.1755-0998.2010.02847.x
Falush D, Stephens M, Pritchard JK (2003) Inference of Population Structure Using Multilocus Genotype Data: Linked Loci and Correlated Allele Frequencies. Genetics 164:1567–1587. https://doi.org/10.1093/genetics/164.4.1567
Fleming WE (1972) Biology of the Japanese beetle, vol 1449. Technical Bulletin. US Department of Agriculture, Washington, DC
Fleming WE (1976) Integrating control of the Japanese beetle: A historical review, vol 1545. US Department of Agriculture, Agricultural Research Service
Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R (1994) DNA primers for amplification of mitochondrial Cytochrome C oxidase subunit I from diverse metazoan invertebrates. Mol Mar Biol Biotechnol 3:294–299
Fox H (1932) The distribution of the Japanese beetle in 1930 and 1931, with special reference to the area of continuous infestation. J Econ Entomol 25:396–407
Fu Y-X (1997) Statistical Tests of Neutrality of Mutations Against Population Growth. Hitchhiking and Background Selection Genetics 147:915–925. https://doi.org/10.1093/genetics/147.2.915
Harry M, Solignac M, Lachaise D (1998) Molecular Evidence for Parallel Evolution of Adaptive Syndromes in Fig-BreedingLissocephala(Drosophilidae). Mol Phylogenet Evol 9:542–551. https://doi.org/10.1006/mpev.1998.0508
Hubisz MJ, Falush D, Stephens M, Pritchard JK (2009) Inferring weak population structure with the assistance of sample group information. Mol Ecol Resour 9:1322–1332. https://doi.org/10.1111/j.1755-0998.2009.02591.x
Hungate BA, Kearns DN, Ogle K, Caron M, Marks JC, Rogg HW (2016) Hydrogen isotopes as a sentinel of biological invasion by the Japanese beetle, Popillia japonica (Newman). PLoS ONE 11:e0149599
Klein M (2008) Popillia japonica (Japanese beetle). Invasive Species Compendium. CAB International, Wallingford, UK. https://doi.org/10.1079/ISC.43599.20203373917
Leigh JW, Bryant D (2015) popart: full-feature software for haplotype network construction. Methods Ecol Evol 6:1110–1116. https://doi.org/10.1111/2041-210X.12410
Marianelli L et al (2019) Long-lasting insecticide-treated nets: A new integrated pest management approach for Popillia japonica (Coleoptera: Scarabaeidae). Integr Environ Assess Manag 15:259–265. https://doi.org/10.1002/ieam.4107
Marianelli L et al (2018) Entomopathogenic nematodes as potential biological control agents of Popillia japonica (Coleoptera, Scarabaeidae) in Piedmont Region (Italy). J Appl Entomol 142:311–318. https://doi.org/10.1111/jen.12470
Mazza G et al (2017) Hexamermis popilliae n. sp.(Nematoda: Mermithidae) parasitizing the Japanese beetle Popillia japonica Newman (Coleoptera: Scarabaeidae) in Italy. Syst Parasitol 94:915–926. https://doi.org/10.1007/s11230-017-9746-0
Mori N, Santoiemma G, Glazer I, Gilioli G, Ciampitti M, Cavagna B, Battisti A (2022) Management of Popillia japonica in container-grown nursery stock in Italy. Phytoparasitica 50:83–89. https://doi.org/10.1007/s12600-021-00948-2
Nardi F, Carapelli A, Dallai R, Roderick GK, Frati F (2005) Population structure and colonization history of the olive fly, Bactrocera oleae (Diptera, Tephritidae). Mol Ecol 14:2729–2738. https://doi.org/10.1111/j.1365-294X.2005.02610.x
Orlov YL, Potapov VN (2004) Complexity: an internet resource for analysis of DNA sequence complexity. Nucleic Acids Res 32:W628–W633. https://doi.org/10.1093/nar/gkh466
Paoli F et al (2017) Differential susceptibility of Popillia japonica 3rd instars to Heterorhabditis bacteriophora (Italian strain) at three different seasons. Biocontrol Sci Technol 27:439–444. https://doi.org/10.1080/09583157.2017.1285866
Pavesi M (2014) Popillia japonica specie aliena invasiva segnalata in Lombardia. L’informatore Agrario 32:53–55
Potter DA, Held DW (2002) Biology and management of the Japanese beetle. Annu Rev Entomol 47:175–205. https://doi.org/10.1146/annurev.ento.47.091201.145153
Pritchard JK, Stephens M, Donnelly P (2000) Inference of Population Structure Using Multilocus Genotype Data. Genetics 155:945–959. https://doi.org/10.1093/genetics/155.2.945
Ramasamy RK, Ramasamy S, Bindroo BB, Naik VG (2014) STRUCTURE PLOT: a program for drawing elegant STRUCTURE bar plots in user friendly interface. SpringerPlus 3:431. https://doi.org/10.1186/2193-1801-3-431
Raymond M, Rousset F, Rousset F (2001) 1995
Regione P (2021) D.D. n. 872 del 12 Ottobre 2021, "Aggiornamento dell'area delimitata per la presenza di Popillia japonica Newman in Piemonte" http://www.regione.piemonte.it/governo/bollettino/abbonati/2021/41/attach/dda1700000872_10400.pdf
Roderick GK, Navajas M (2003) Genes in new environments: genetics and evolution in biological control. Nat Rev Genet 4:889–899. https://doi.org/10.1038/nrg1201
Rousset F (2008) genepop’007: a complete re-implementation of the genepop software for Windows and Linux. Mol Ecol Resour 8:103–106. https://doi.org/10.1111/j.1471-8286.2007.01931.x
Scott MJ et al (2020) Genomic analyses of a livestock pest, the New World screwworm, find potential targets for genetic control programs. Commun Biology 3:424. https://doi.org/10.1038/s42003-020-01152-4
Simon C, Frati F, Beckenbach A, Crespi B, Liu H, Flook P (1994) Evolution, Weighting, and Phylogenetic Utility of Mitochondrial Gene Sequences and a Compilation of Conserved Polymerase Chain Reaction Primers. Ann Entomol Soc Am 87:651–701. http://doi.org/10.1093/aesa/87.6.651
Tajima F (1989) Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123:585–595. https://doi.org/10.1093/genetics/123.3.585
Torrini G et al (2020) Evaluation of indigenous entomopathogenic nematodes as potential biocontrol agents against Popillia japonica (Coleoptera: Scarabaeidae) in Northern Italy. Insects 11:804. https://doi.org/10.3390/insects11110804
Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, Remm M, Rozen SG (2012) Primer3—new capabilities and interfaces. Nucleic Acids Res 40:e115–e115. https://doi.org/10.1093/nar/gks596
Ye J, Coulouris G, Zaretskaya I, Cutcutache I, Rozen S, Madden TL (2012) Primer-BLAST: A tool to design target-specific primers for polymerase chain reaction. BMC Bioinformatics 13:134. https://doi.org/10.1186/1471-2105-13-134

Table 1 Primers used in this work and annealing temperature used. Microsatellite primer pairs were annotated with parameters of their reference allele.

Locus	Locus Type	Primer Forward	Forward Primer Sequence	Primer Reverse	Reverse Primer Sequence	Tm	Reference repeat	Length range	Alleles	Missing data	He
A	Microsatellite	Pj_35808_Fw	TGTTTGTCCTCGGCGATGC	Pj_35808_Rv	ACCACCGTCTAACAGCAGG	56° C	(GT30)	142-232	29	0%	0.87
B	Microsatellite	Pj_14476_Fw	ACAGAAGGTCCGATAGCATGG	Pj_14476_Rv	CCATGCCCTATCAAAGTCAAGC	56° C	(AG)13	232-268	11	0%	0.69
C	Microsatellite	Pj_3266_Fw	GCTGCAACTCTGTTCCATTTCC	Pj_3266_Rv	TAGCACCTTACTCTCCAGGC	56° C	(TA)17	228-246	10	0.9%	0.58
D	Microsatellite	Pj_6327_Fw	CGACAACGCAAGCAAGTCC	Pj_6327_Rv	AGTAGCATGGAATGACTGACGG	56° C	(TAA)15	236-275	6	4.5%	0.34
E	Microsatellite	Pj_6701_Fw	ACAGTGCTTTTCGTTACGATGC	Pj_6701_Rv	CGTTTGAAGCCGTGATAGGG	56° C	(AAT)14	230-287	17	8.9%	0.84
G	Microsatellite	Pj_445_Fw	CCGATCCAATCGTAATGGAAGC	Pj_445_Rv	CTTTGAAGTAGCGGTAGTAGGC	56° C	(TG)15	156-178	10	17.9%	0.62
H	Microsatellite	Pj_24134_Fw	TAAACCCACCGGAAGCTTGG	Pj_24134_Rv	AGTGAGTGCCCTCGTAGACC	56° C	(TC)14	166-200	17	23.7%	0.79
I	Microsatellite	Pj_31638_Fw	CGTGTCCATTCGGATCGTTAG	Pj_31638_Rv	ACAGCCTTGGAAAATGGACAG	54° C	(AAT)22	153-249	27	42.0%	0.91
L	Microsatellite	Pj_20327_Fw	GCAATTGGTAGAGAACCTCTGG	Pj_20327_Rv	GAAGCCTGCATTGGAGGAAG	56° C	(TAA)26	153-240	25	20.5%	0.85
COXI	Mitochondrial	LCO1490	GCTCAACAAATCATAAAGATATTGG	HCO2198	TAAACTTCAGGGTGACCAAAAAATCA	45° C	-	-	-	-	-
CytB	Mitochondrial	CP1	GATGATGAAATTTTGGATC	CB-N-11367	ATAACTCCTCCTAATTTATTAGGAAT	40° C	-	-	-	-	-

Table 2 Genetic diversity indices and Fu’s Fs and Tajima’s D tests in COX I and CytB gene fragments of Popillia japonica. Abbreviations: n – Unique haplotype number; Hd - haplotype diversity, S - Number of polymorphic sites; η(S) - Total number of singleton mutations; k - average number of pairwise nucleotide differences, π - nucleotide diversity per site

Number of positions

η(S)

Fu’s F_S

(p ≤ 0. 05)

Tajima’s D

(p ≤ 0.05)

COX I

657

0.828 ± 0.020

5.35477

0.00815 ± 0.00100

-8.23268

p = 0.042

-1.14028

p = 0.099

CytB

744

0.870 ± 0.019

7.04582

0.00947 ± 0.00082

-18.89330

p = 0.001

-1.42761

p = 0.050

S1.xlsx
Table S1 List of all samples used in this work with sampling dates and sites, geographical coordinates, mitochondrial haplotype found and sample groups assigned for microsatellite analysis
S2.pdf
Fig S2 Microsatellite analysis of sub-datasets. a: Fig 1a repeated to allow comparisons; b: Japan only, 2 groups; c: Japan and USA, 2 groups; d: Japan and USA, 4 groups; e: USA, Azores and Italy, 2 groups; f: USA, Azores and Italy, 4 groups.
S3.tif
Fig S3 Single tree from COX I obtained in Bayesian and Maximum likelihood analyses. Posterior Probabilities and Bootstrap values are labeled above the appropriate internodes.
S4.tif
Fig S4 Single tree from CytB obtained in Bayesian and Maximum likelihood analyses. Posterior Probabilities and Bootstrap values are labeled above the appropriate internodes.
SupplementaryData.str

Download PDF

Version 1

posted

You are reading this latest preprint version

Tracing the dispersal route of the invasive Japanese beetle Popillia japonica

Status:

Version 1

Abstract

Figures

Key Messages

Introduction

Materials And Methods

Sampling and DNA extraction

Microsatellite development

Optimization of candidate loci

Genotyping

Microsatellite analysis

Amplification of mitochondrial targets, sequencing and data analysis

Results

Dataset evaluation and basic statistics

Microsatellite data

Mitochondrial data

Discussion

Conclusions

Declarations

References

Tables

Supplementary Files

Status:

Version 1