Microsatellite and Morphological Analyses Reveal Unexpected Diversity in Lymantria dispar in China

Abstract

The gyspy moth Lymantria dispar Linnaeus, a widely distributed leaf-eating pest, is considered geographically isolated in the world, with two Asian gypsy moth subspecies, Lymantria dispar asiatica and Lymantria dispar japonica. In China, only one subspecies, L. d. asiatica, has been observed. In this study, we characterized gypsy moth diversity and divergence using 427 samples covering a wide range of the species distribution, with a focus on sampling along a latitudinal gradient in China. We combine the quantitative analysis of male genitalia and the genetic diversity analysis of nine microsatellite loci of nuclear genes nuclear genes to study the structure of gypsy moth individuals in 23 locations in the world and the male genitalia of gypsy moths in some areas. In mixed ancestry model-based clustering analyses based on nuclear simple sequence repeats, gypsy moths were divided into three well-known subspecies, a unique North American cluster, and a southern Chinese cluster with differentiation between the Asian gypsy moth and European gypsy moth. We also found individuals identified as European gypsy moths in two distant regions in China. The results of a quantitative analysis of male genitalia characteristics were consistent with an analysis of genetic structure and revealed the differentiation of gypsy moths in southern China and of hybrids suspected to be associated with L. d. japonica in the Russian Far East. Admixture in gypsy moths can be explained by many factors such as human transport. In China, we detected European gypsy moths, and found unexpectedly high genetic diversity within populations across a wide range of latitudes.

1. Introduction

Genetic lineages within species often exhibit striking geographic patterns, and the genetic structure of populations is shaped by various processes [1,2,3]. Many species inhabit a wide range of latitudes. Accordingly, integrative analyses are important for the characterization of diversity and species discovery [4,5,6,7]. In most studies, the difference in morphology, such as external morphological characters and genitalia, is the first and most noticeable part of the research interest. In addition, the development of molecular technology has brought new ideas to the study of species. Increasing evidence demonstrates the importance of considering both morphological and genetic variation in studies of population structure and species differentiation [8,9,10,11,12,13,14].

Lymantria dispar Linnaeus (Lepidoptera: Erebidae: Lymantriinae), also known as the gypsy moth, is a widely distributed omnivorous forest pest. The gypsy moth is native to Eurasia, with a wide distribution spanning the Holarctic Sweden, Norway, the Russian Far East, Africa, Japan, and North America [15,16]. It has been recorded in Hebei, Shanxi, Inner Mongolia, Liaoning, Jilin, Heilongjiang, Shandong, Hubei, Hunan, Shaanxi, Gansu, and Xinjiang in China [17]. L. dispar is not an indigenous species in North America; from 1868 to 1869, it was artificially imported from France and escaped in Massachusetts, expanding rapidly to Canada [18,19]. L. dispar was divided into three typical subspecies by morphological characteristics, female flight ability, and geographical region [20]: L. dispar dispar in Europe and North America, also known as the European gypsy moth (EGM); and L. dispar asiatica and L. dispar japonica, collectively referred to as the Asian gypsy moth (AGM). The United States Department of Agriculture (USDA) believes that female AGM have greater flight ability than males and are able to fly up to 25 miles (40 kilometers) in some cases. The wide range of host plants and the strong flight ability of females allow AGM to rapidly enter uninfected areas and cause damage. Since 1991, the North American Plant Protection Organization (NAPPO) and its member states, as well as Chile and New Zealand, quarantine vessels from AGM habitats at specific times [21,22].

Owing to the damage caused by L. dispar to forest trees and its wide distribution, the population structure of the species based on molecular markers has been a focus of research. Since mitochondrial DNA (mtDNA) revealed the genetic inheritance of the NAGM and gypsy moth in Europe; and haplotypes differ between Asian and North American groups [23]. To quickly and accurately distinguish between AGMs and EGMs, many molecular marker techniques have been explored, widely used in quarantine departments [24,25,26]. It was proposed that gypsy moths spread from east to west in Eurasia [27,28]. An analysis of mitochondrial and nuclear genes of worldwide gypsy moths divided gypsy moths into three groups: North America, Europe/Siberia, and Asia. Gypsy moths could be divided into four clusters corresponding to the three subspecies and populations in North America [29,30]. The Korean Peninsula and its adjacent areas have mixed ancestry, with populations that are not concordant with genetic differentiation among subspecies [31]. A recent analysis showed that the gypsy moth could be divided into three lineages: a Northeast Asia and Japan lineage; a Europe and Central Asia lineage; and a Trans-Caucasus lineage, potentially explaining the diffusion of the species. Lymantria umbrosa (Butler) is a sister group to other gypsy moth populations [32]. The discordance between nuclear and mitochondrial gene analyses should be resolved in future studies, due to the different evolution rate [33].

In these studies, gypsy moths in China are classified as the typical AGM (L. dispar asiatica) based on single sampling locations, mostly in northern China. However, genetic differentiation between northern and southern China has been neglected, despite complex geographical and climatic conditions. Some recent studies accounted for this potential complexity of lineages in China. An analysis of different fragments of the mitochondrial gene in different regions spanning the majority of the L. dispar habitat in China identified a few locations in China with haplotypes that are closely related to EGM haplotypes. Some gypsy moth samples from southern China form a sister group with typical AGMs in a phylogenetic tree, and the genetic structure of the gypsy moth is complex [34,35].

Although mitochondrial genes have various advantages for population genetic analyses, including maternal inheritance, easy amplification, and low molecular weight. The rate of evolution of nuclear genes is slower than that of mitochondrial genes, but these genes are more stable and more closely related to the regulation of biological traits [36,37]. Of note, there is the potential for mito-nuclear discordance when assessing gypsy moth population genetic structure. In particular, in some populations of L. dispar dispar from Central Asia based on mitochondrial haplotypes, females often display AGM-like flight ability, and their nuclear genomes appear to have more in common with those of L. dispar asiatica than L. dispar dispar [29,30,38,39,40].

In previous studies, gypsy moth sampling in China has been insufficient, and analyses have been limited to mitochondrial genes. Despite high mitochondrial gene diversity in the Chinese gypsy moth, similar trends in nuclear genes and morphology have not been demonstrated. To better probe genetic differentiation in gypsy moths in China, we used more comprehensive sample collection, including collection of subspecies found in other regions (L. d. dispar and L. d. japonica), combined with nuclear genetic analyses and morphological evaluations of male genitalia, to explore gypsy moth diversity.

2. Materials and Methods

2.1. Specimen Collection and DNA Extraction

Adult male gypsy moths were collected by pheromone traps. Additionally, egg masses collected from the field were raised to adults by artificial diet feeding technology in the laboratory [41]. Two widely recognized gypsy moth groups were included: AGMs and EGMs. In total, 427 individuals were sampled from 23 regions in seven countries (Figure 1).

The samples included three typical gypsy moth subspecies, L. dispar dispar, L. dispar asiatica, and L. dispar japonica. More detailed information for all samples is provided in

Table 1. Sample information for Lymantria dispar.

Number	Code	Sampling Location	Longitude	Latitude	N1	N2
1	HLJ	Hegang, Heilongjiang †	130°49′44″E	47°34′39″N	20	10
2	LN	Wafangdian, Liaoning	121°58′44″E	39°37′37″N	10
3	BJ	Beijing	116°12′25″E	40°18′01″N	10
4	HB	Zhangjiakou, Hebei †	115°52′23″E	41°00′11″N	20	10
5	CRS	Charisuzhen, Tongliao, Inner Mongolia	123°28′34″E	43°07′36″N	20
6	Kdr	Hulunbeir, Inner Mongolia †	120°42′43″E	49°17′08″N	20	10
7	TL	Tongliao, Inner Mongolia	122°15′47″E	43°37′02″N	20
8	XJ	Urumqi, Xinjiang	87°56′48″E	44°10′06″N	17
9	AH	Lu’an, Anhui †	116°30′27″E	31°45′10″N	20	10
10	SX	Yuncheng, Shanxi	111°00′14″E	35°01′21″N	10
11	CD	Chengdu, Sichuan †	104°03′58″E	30°34′23″N	20	10
12	YN	Lushui, Yunnan †	98°51′15″E	25°51′03″N	20	10
13	XF	Xifeng, Guizhou †	106°44′25″E	27°05′25″N	20	10
14	ZY	Zunyi, Guizhou †	106°56′14″E	27°42′23″N	20	4
15	RM	Primorsky Krai, Russia †	132°00′40″E	43°11′11″N	20	10
16	Hu	Honshu, Japan †	139°02′35″E	35°04′28″N	20	10
17	RBI	Shira, Russia †	89°57′52″E	54°30′06″N	20	10
18	MA	Massachusetts, USA	71°50′02″W	42°30′30″N	20
19	NC	North Carolina, USA	83°14′36″W	35°42′06″N	20
20	CT	Connecticut, USA †	72°41′27″W	41°37′39″N	20	10
21	JL	Lithuania	25°14′12″E	54°41′09″N	20
22	KG	Greece †	22°41′10″E	39°48′06″N	20	10
23	FR	France †	02°43′58″E	48°52′06″N	20	10

1–14: Localities in China. N₁: Number of samples in the microsatellite analysis, N₂: Number of samples in the morphological analysis (locations of samples used in morphological analyses are marked with †).

.

After drying all adult samples, about 20 mg thoracic muscle tissue was removed and ground into a powder in liquid nitrogen. Genomic DNA was extracted and purified using the Rapid and Non-toxic Insect Genome DNA Extraction Kit (Demeter DNA, Beijing, China). The quality of extracted DNA was tested using the NanoDrop 8000 (Thermo Fisher Scientific, Waltham, MA, USA). When the 260/280 ratio reached ~1.8, the DNA quality was considered good [42] and further experiments were performed. All dried samples and extracted DNA were stored at −20 °C.

2.2. Microsatellite Loci and PCR Amplification

Information for many microsatellite loci was obtained from previous studies [28,30,43]. A small number of DNA samples from different regions were used to screen all the primers we found. The DNA fragments based on nine microsatellite primers were successfully amplified, and these microsatellite primers (

Table 2. Summary of microsatellite loci.

Locus Name	Motif	Label	Forward Primer	Reverse Primer	Annealing Temperature (°C)
10F1	(AC)	FAM	CGCACAAAGCTCTCAGATGA	CGTTACCGCGTGTCTAGATT	53
49	(TGA)	HEX	GAAGCCTACATTCAGCAGTTG	GAAATCCGTCCATCCATTTG	54
101	(TGA)	ROX	AATTTACCCTTGCGTTATGTAGAC	ACATATTCGAACAGTTGTTTCATAA	54
106	(TGA)	ROX	AGGCTCGATGCCAGTAGTGG	ACAAAGCCAATCGGATAGAACA	54
107	(GT)	ROX	TCTGAAGCGAGATGAACTGG	TAAGCTTAGACCTCCTCCAG	53
138	(GT)	HEX	TTCGTTCAGTGAGCGAGAGA	CTCCATACCCCAATCAAGAC	54
202	(AC)	FAM	TCCCATATCTGTCCACACCA	AATCCATTAAAATCGGTCTAGCC	53
238	(CA)	HEX	ACTGTTCGTTTATTCAATAGTGTTGG	ATATCCCTTAGTCGCCTTTTACG	54
254	(CA)	ROX	TACTGTTTGAAGTCGGTTTTGC	GATGACTAGCGTATTCAATACGCA	54

Primer source: 10F1, 49, and 202 from Bogdanowicz et al. [28]; 107 and 128 from Koshio et al. [43]; 101, 106, 238, and 254 from Wu et al. [30].

, Supplementary Figure S1) were ultimately included in the study. The overall Ewens Watterson test for neutrality was carried out on nine microsatellite loci. The results showed that all the loci adopt the neutral test, and there was no case that loci could not be used because of deviating from the neutral test (Supplementary Table S1).

Polymerase chain reaction was performed using a 20 μL volume containing 7 μL deionized water, 10 μL of 2× T5 Super PCR Mix (PAGE) (Tsingke, Beijing, China), 1 μL forward primer, 1 μL reverse primer, and 1 μL DNA template (30–60 ng/μL). Fluorescent labels, including ROX (red), HEX (green), and FAM (blue) (Tsingke), were used to label the forward primers in each pair. The ROX-labeled fragment had the shortest length and the FAM-labeled fragment had the longest; this length difference prevented errors. The PCR cycling conditions for microsatellite loci 49, 101, 106, 138, 238, and 254 included pre-denaturation at 98 °C for 2 min, followed by denaturation at 95 °C for 10 s, renaturation at 54 °C for 10 s, and extension at 72 °C for 15 s for 32 cycles, with a final extension at 72 °C for 15 min. The cycling conditions for 10F1, 107, and 202 were the same, but the annealing temperature was 53 °C. These procedures were performed in dark conditions to avoid quenching of the fluorescent labels. The PCR products were stored at 4 °C, and sequencing was performed by Tsingke Biological Technology Co., Ltd. using the ABI3730xl capillary electrophoresis sequencer, with GS500 as an internal standard.

Fragment length analyses and genotyping were performed using GeneMarker v.2.2.0 [44]. Samples that failed to successfully amplify all 9 primers and had inappropriate DNA concentrations were excluded from this experiment, and the number of samples of each population for microsatellite analysis was shown in Table 1. CONVERT v.1.3.1 was used to convert diploid genotype data files into the formats that can be directly read by most programs for population genetic analyses [45]. Micro-Checker v.2.2.3 was used to identify genotyping errors caused by null alleles, short allele dominance (large allele dropout), and scoring errors due to stuttering [46].

2.3. Analysis of Genetic Diversity

Genepop v.4.7.0 [47] was used to evaluate deviations from Hardy–Weinberg equilibrium (HWE) and linkage disequilibrium (LD) by the Markov chain method, with 10,000 dememorizations and 50 million iterations. Sequential Bonferroni-corrected p-values (p < 0.0153) were used for each population and locus in HWE equilibrium. FSTAT 2.9.4 [48] was used to calculate the average genetic diversity (D), number of alleles (N), average number of alleles, mean number of alleles per locus per population (MNa), inbreeding coefficient (F_IS), mean allelic richness after correcting for sample size differences (Rs), and allelic richness using a minimum sample size of ten. The Excel add-in MS_tools [49] was used to calculate the observed heterozygosity (H_O), expected heterozygosity (H_E), and polymorphism information content (PIC) [50]. PIC < 0.25 indicated low polymorphism; 0.25 < PIC < 0.5 indicated an intermediate level of polymorphism; PIC > 0.5 indicated high polymorphism. Genetic differentiation, as estimated by the fixation index (F_ST), among populations was calculated using Genetix v.4.05 [51], and significance was evaluated by 1000 bootstrap replicates. GenALEX 6.5 [52,53] was used to calculate the number of private alleles over all loci (P_A) and the pairwise population matrix of Nei’s unbiased genetic distances. Isolation-by-distance was evaluated using GenALEX 6.5 by the Mantel test [54,55] of the correlation between the geographic distance matrix and Nei’s unbiased genetic distance matrix.

2.4. Population Genetic Structure

The genetic structure of 23 populations of L. dispar was inferred from microsatellite data using admixture models applied in STRUCTURE v.2.3 [56]. A genetic clustering analysis was performed using 427 gypsy moth individuals and nine loci. The admixture model was selected to allow mixed pedigrees. K-values from 1 to 23 were considered, the burn-in period was set to 50,000, and 200,000 formal Markov chain Monte Carlo simulations (MCMC) were run, with 20 iterations for each K-value.

To select the most likely k-value [57], the STRUCTURE results were evaluated using Structure Harvest [58], and the appropriate grouping was selected by considering the Ln Pr(X|K) value and the ΔK value. When the maximum ΔK value could not fully express the difference between populations, the grouping result for the next ΔK value was considered. After determining the grouping results, results for 20 iterations using the optimal k-value were averaged using CLUMPP v.1.1.2 [59], and the stability of the optimal k-value was evaluated by H’, using 0.9 (indicating that group similarity exceeded 90%) as a threshold, and average results were mapped. To fully characterize the complex genetic structure, the individual membership coefficients < 0.05 were not simplified.

Gypsy moth individuals were grouped by the K-value calculated using STRUCTURE v.2.3, and genetic variance between populations, individuals, and groups was evaluated by analysis of molecular variance (AMOVA) implemented in ARLEQUIN v.3.5 [60]. The significance of the covariance at each level was calculated using 5000 permutations.

Bottleneck v.1.2.02 was used to detect bottleneck effects in 23 populations of L. dispar. The calculations were performed using 1000 iterations with three statistical tests (i.e., the Wilcoxon signed rank test, sign test, and standardized difference test) and with three mutation models (Infinite Allele Model, IAM; Stepwise Mutation Model, SMM; Two-phased model, TPM). The Wilcoxon signed rank test is preferable to the sign test or the standardized difference test because it requires fewer sites and has higher power; for 15–40 samples and 10–15 polymorphic loci, credible test results can be obtained [61,62,63]. Additionally, most microsatellite data are more effectively analyzed using the TPM model than the IAM and SMM models [63]. Therefore, in the analysis of bottlenecks for nine polymorphic loci in 23 populations of gypsy moth, the TPM model results using the Wilcoxon signed rank test are likely the most realistic.

2.5. Measurements of Morphological Characteristics and Statistical Analyses

For analyses of morphological characteristics, based on the genetic structure of the genes, we selected gypsy moths from the habitats of the three subspecies (L.d.dispar, L.d.asiatica and L. d. japonica) and southern China and North America. Nine characters related to the genitalia of male gypsy moths were compared among 14 geographic populations (Figure 2). A total of 134 male samples were dissected (Table 1). Multiple samples were randomly sampled in the population. Samples used for analyses of morphological characteristics were basically obtained from the samples used in molecular experiments.

Our approach for anatomical analyses of gypsy moth genitalia was based on these previous studies [10,64,65]. Moreover, according to the characteristics of insects, some improvements have been made. The final 5 mm of the abdomen was cut and placed in a tube with 10% KOH solution for 1 to 2 h in a water bath to dissolve impurities outside of the genitalia. When the abdominal section was visibly translucent, it was removed and immersed in a glass dish containing distilled water. Under a MSZ745 stereo-microscope (Nikon, Shanghai, China) the somite was dissected from the side using fine tweezers to remove excess undissolved fat and other tissues. The genitalia were removed and placed in a concave porcelain plate with Bouin’s fixative solution for 2 to 4 h. Because the fixative has a staining effect and its own color, gypsy moth genitalia are yellow to brownish-yellow, genitalia could be clearly observed under a microscope, and no additional dyes were used. The fixed material was rinsed in 90% alcohol, 75% alcohol for 30 min, and 75% alcohol, followed by storage in a small centrifuge tube at 4 °C.

The treated genitalia were placed on a glass slide, the valva was stretched, the posture was adjusted, and another glass slide was used to gently flatten the specimen. A M205FA Fluorescence Automated Microscope (Leica, Wetzlar, Germany) was used to obtain images of specimens, and characters were measured using Leica Application Suite (LAS). Statistical analyses were performed using RStudio [66]. A principal component analysis (PCA) was used to summarize variation in male genitalia. A correspondence analysis (CA) was used explore the relationships between characters and specimens. Combining the PCA loadings and CA results, characters with a substantial influence on the specimen were identified and used for a new PCA to reduce noise.

3. Results

3.1. Population Genetic Diversity and Differentiation

Genetic diversity for 427 gypsy moth samples from 23 regions based on nine microsatellite loci are summarized in

Table 3. Genetic diversity in 23 populations of Lymantria dispar.

CODE	Na	AP	Rs	D	mPIC	HE	HO	FIS	HWE	TPM
HLJ	7.33	3	5.904	0.755	0.695	0.732	0.572	0.243	0.0000 *	0.8203
LN	5.56	2	5.556	0.714	0.636	0.672	0.589	0.175	0.0016 *	1.0000
BJ	4.89	2	4.889	0.644	0.570	0.614	0.678	−0.052	0.2384	0.7344
HB	6.44	1	5.216	0.666	0.609	0.647	0.561	0.157	0.0001 *	0.8203
CRS	6.67	2	5.304	0.654	0.601	0.635	0.550	0.158	0.0005 *	0.6523
Kdr	4.44	0	3.943	0.638	0.556	0.619	0.550	0.137	0.0200	0.0137 ‡
TL	7.33	3	5.867	0.709	0.659	0.688	0.578	0.185	0.0000 *	0.7344
XJ	4.89	2	4.493	0.686	0.601	0.661	0.516	0.247	0.0000 *	0.0098 ‡‡
AH	7.22	3	5.749	0.700	0.649	0.680	0.611	0.127	0.0002 *	0.6523
SX	5.56	0	5.556	0.712	0.631	0.672	0.633	0.110	0.0255	0.3008
CD	7.22	5	5.863	0.746	0.688	0.724	0.589	0.211	0.0000 *	0.5703
YN	3.78	4	3.452	0.523	0.452	0.508	0.439	0.162	0.0020 *	0.4258
XF	5.44	0	4.850	0.720	0.647	0.697	0.544	0.243	0.0000 *	0.0195 ‡
ZY	6.67	0	5.568	0.710	0.652	0.689	0.567	0.202	0.0000 *	0.1641
RM	3.67	1	3.372	0.539	0.470	0.526	0.533	0.011	0.2806	0.2500
Hu	2.89	0	2.646	0.493	0.403	0.479	0.439	0.110	0.0679	0.1641
RBI	3.33	2	3.072	0.489	0.426	0.475	0.444	0.090	0.1005	0.5469
MA	2.78	0	2.559	0.479	0.382	0.465	0.422	0.118	0.0318	0.1289
NC	2.89	0	2.517	0.469	0.368	0.453	0.283	0.396	0.0000 *	0.2031
CT	2.89	0	2.708	0.488	0.404	0.473	0.378	0.226	0.0002 *	0.0039 ‡‡
JL	3.33	0	3.064	0.501	0.429	0.487	0.422	0.158	0.0733	0.3594
KG	3.00	0	2.579	0.367	0.311	0.356	0.311	0.152	0.0491	0.6406
FR	5.22	0	4.564	0.649	0.587	0.631	0.600	0.075	0.0018 *	1.0000

Na: mean number of alleles per locality; A_P: number of private alleles over all loci; Rs: allelic richness per population based on minimum sample size of ten diploid individuals; D: gene diversity per population; mPIC: mean polymorphism information content per population over all loci; H_E: expected heterozygosity; H_O: observed heterozygosity; F_IS: inbreeding coefficient; HWE: Hardy–Weinberg equilibrium (* statistical significance for departure from HWE after sequential Bonferroni correction); TPM (^‡ indicates p < 0.05 by a two-tailed model).

.

Using Micro-Checker, we found that only a few loci showed possible null alleles in a few populations, and the frequency of null alleles was low. A few loci in some populations showed significant deviation from HWE (H_O < H_E), and samples from 14 localities showed deviations from HWE in a comprehensive analysis (Supplementary Table S2). Localities with samples exhibiting relatedness were similar to those with samples deviating from HWE. We did not detect relatedness shared among all populations or a large number of loci in a population. Based on these analyses, we believe that the observed deviations from HWE could be attributed to relatedness and null alleles.

We detected 128 alleles for nine microsatellite loci, with an average of 14.2 alleles per loci (range, 7–20 alleles). The mean number of alleles at different locations ranged from 2.78 to 7.33. The mean number of alleles was highest in the HLJ and TL populations in northern China and lowest in the MA population. The number of alleles in gypsy moth populations in China was generally higher than that in foreign populations. Private alleles were mainly detected in samples from China and Russia in East Asia, and the number of private alleles ranged from one to five. Samples collected in Chengdu in southern China had the most private alleles, followed by gypsy moths in Yunnan, with four private alleles. After correcting for differences in sample size, the allelic richness in gypsy moths from different localities ranged from 2.517 to 5.904, with a mean value of 4.317. Allelic richness was below average in Europe, North America, and Japan, except in France. Gene diversity is also part of total genetic diversity, and sampling has a relatively low impact on this parameter. Compared with F-statistics, focused on the effects of single loci in the total population and subpopulations, gene diversity is more concerned with genetic variation with many loci and the total genome [67]. As shown in Table 3, the mean gene diversity was 0.611. Gene diversity was highest in Heilongjiang and lowest in Greece. For nine loci in each population, the average PIC was 0.311 to 0.695, exceeding 0.25. High polymorphism levels were observed in different populations. The PICs for nine loci were higher than the average value in HLJ, LN, TL, AH, SX, CD, XF, and ZY (mPIC > 0.7). Since expected heterozygosity is calculated assuming HWE under ideal conditions, it is generally higher than observed heterozygosity [68]. The observed heterozygosity of all gypsy moths was 0.591, indicating high genetic diversity.

The inbreeding coefficient F_IS refers to the deviation between observed and expected heterozygosity within a population and is an indicator of the degree of inbreeding. F_IS = 0 indicates random mating. Positive values indicate some degree of inbreeding, and negative values indicate outcrossing [69,70]. In the 23 gypsy moth populations, F_IS ranged from −0.052 to 0.369 (mean, 0.158), with the highest value observed in North Carolina and the lowest value in Beijing. Inbreeding and heterozygote depletion were detected in the sample set.

When F_ST < 0.05, there is negligible genetic differentiation between populations; when 0.05 < F_ST < 0.15, there is moderate genetic differentiation among populations; when 0.15 < F_ST < 0.25, there is high genetic differentiation among populations; and when F_ST > 0.25, populations are extremely differentiated [71]. As shown in a heat map in Figure 3, pairwise F_ST values ranged from −0.00257 to 0.48286. The F_ST values for comparisons between TL and CRS populations and between HLJ and LN populations were low, indicating a lack of differentiation, and values were highest for the comparison between KG and RM. The degree of genetic differentiation between Chinese populations and other populations was relatively high, and differentiation was low among populations within China. However, the paired F_ST values between Yunnan and other populations were unique. This population exhibited high differentiation from all other populations (including other populations in China). The pairwise F_ST value for YN and KG was 0.44574, indicating extremely high differentiation. The F_ST values for Chinese Xifeng and Kuduer populations, and for European populations indicated that the two populations in China were poorly differentiated from gypsy moths in France and Lithuania.

To explore the influence of isolation-by-distance on phylogenetic relationships, an unbiased genetic distance matrix was compared with a log-transformed geographic distance (GGD) matrix (Figure 4). We detected a positive correlation (R = 0.501, p < 0.05, y = 0.253x + 0.0684, Mantel test), suggesting that geographic distance is related to genetic distance for 23 populations of gypsy moths based on SSRs. Geographic isolation might be a major determinant of genetic differentiation among gypsy moths.

3.2. Genetic Structure and Variation

We performed a cluster analysis with an admixture model using STRUCTURE. We found that when K = 2, the peak value of ∆K was the highest, and the division of populations was the most reasonable. In particular, 23 populations were divided into two clusters (AGM and EGM clusters, represented by pink and blue-violet in Figure 5, respectively). Most populations in China, as well as Honshu Island and the Russian Far East, were assigned to the same group. The Xifeng and Kuduer populations clustered with the European, North American, and Russian Siberian populations. The gypsy moths from China, the Russian Far East, and Japan were assigned to the AGM cluster and were divided according to geographical location. Gypsy moths from Europe and North America were assigned to the EGM cluster.

When K = 2, the samples were effectively divided into two types, but genetic diversity in China was not clearly displayed. Therefore, the second highest peak for ∆K was selected to further explore the genetic structure of the gypsy moth, K = 5. The 23 populations of gypsy moth were divided into five clusters indicated by green, yellow, blue, red, and purple in Figure 6. The North America and Russian Shira populations were assigned to the same group. European specimens and some samples from Xifeng and Kuduer were grouped together (European cluster). A few samples from Xifeng and other populations in southern China were clustered (South China cluster). Although, samples from AH were grouped in South China cluster, but there are about seven samples showing a closer relationship with the northern Chinese cluster. The gypsy moths in northern China and the Primorsky Krai in Russia were assigned to the same group (North China and the Russian Far East cluster). The population in Honshu Island, Japan was not grouped with other populations (Japan cluster).

We divided gypsy moth populations into two groups and five groups according to an admixture model, and used AMOVA to analyze groups and levels (

Table 4. Analysis of molecular variance.

	Source of Variation	d.f.	Sum of Squares	Variance Components	Percentage of Variation (%)
K = 2	Among groups	1	140.066	0.28129 Va	8.01 *
	Among populations within groups	21	461.404	0.52072 Vb	14.83 *
K = 5	Among groups	4	310.208	0.37513 Va	10.87 *
	Among populations within groups	18	291.261	0.36811 Vb	10.66 *
	Within populations	831	2251.018	2.70881 Vc	78.47 *
	Total	853	2852.487	3.45204	— —

* p < 0.01; d.f.: Degrees of freedom.

). When the gypsy moth populations were divided into two groups, genetic variation within groups was 8.01% (p < 0.01), and variation among populations within groups was 14.83% (p < 0.01). When populations were divided into five groups, genetic variation was 10.87% (p < 0.01) within groups and 10.66% (p < 0.01) among populations within groups. Compared with other source of variation, Genetic variation within populations was highest, 78.47% (p < 0.01), but there was also clear divergence between groups.

Under the two-tailed Wilcoxon signed rank test with the TPM model (Table 3), Kuduer (Kdr), Xinjiang (XJ), Xifeng (XF), and Connecticut populations (CT) experienced a bottleneck (p < 0.05), and these bottlenecks were highly significant in the Xinjiang and Connecticut populations (p < 0.01).

3.3. Morphological Analysis

Pogue described subspecies morphologies of L. dispar in detail in 2007 and showed that male genitalia of the three subspecies are significantly different with respect to the ratio of the male phallus to the genital capsule [20]. Therefore, after identification according to the morphological description, we performed a detailed analysis of the morphology of male genitalia of the gypsy moth among populations.

We selected nine characters of male genitalia (LUNC, length of the uncus; LSAC, length of the sacculus; LPHA, length of the phallus; WUNC, width of the uncus; LVAL, length of the valva; WVAL, width of the valva; WJUX, width of the juxta; LUTS, length from the uncus to sacculus; WSAC, width of the sacculus) for a quantitative morphological analysis. The results of a PCA of these variables are summarized in

Table 5. Principal component analysis of male genitalia of L. dispar.

	Component 1a	Component 2a	Communality	Component 1b	Component 2b
Eigenvalue	5.723	1.209		3.292	1.028
Standard deviation	2.392	1.100		1.814	1.013
Proportion of variance	0.636	0.134		0.658	0.205
Cumulative proportion	0.636	0.770		0.658	0.864
Loadings:
LUNC	0.371	0.109	0.802	0.503	0.150
LSAC	0.352	0.109	0.722
LPHA	0.347	0.413	0.895	0.465	0.396
WUNC	0.219	−0.504	0.581
LVAL	0.387	0.162	0.614
WVAL	0.320	−0.150	0.889	0.444	−0.201
WJUX	0.351	−0.114	0.722
LUTS	0.407	0.000	0.955	0.539
WSAC	0.170	−0.696	0.751	0.208	−0.883

Component 1a, Component 2a: Compositions of principal component analysis of 9 characters; Component 1b, Component 2b: Components of the principal component analysis of the four characters, after noise reduction.

.

The first two principal components accounted for 63.6% and 13.4% of the variance, and the cumulative contribution rate reached 77%. Four characteristics (LUNC, LPHA, WVAL, and LUTS) had the highest load (p > 80%).

In a PCA of 14 populations of gypsy moth with nine complete characters, we found little differentiation, except for samples from Japan, with relatively wide scatter for points in the same area (Figure 7A). Therefore, we performed a CA of the relationship between samples and features (Figure 7B). After reducing the characteristics, male genitalia of samples from Japan and the Russian Far East were closely related, and the genitalia were relatively large, with most differentiation in the Japanese population. Genitalia in the YN, CT, and ZY populations were smaller, and samples in other regions were basically indistinguishable under the two principal components.

4. Discussion

4.1. Genetic Variation in Gypsy Moth Populations

Genetic diversity plays a vital role in species evolution. Without considering other factors, genetic diversity increases the population growth rate, and the degree of heterozygosity is related to fitness. Bottleneck effects and inbreeding could alter genetic diversity and affect species persistence under stressful conditions [72,73,74]. The geographical division among the three major subspecies of L. dispar is unclear. The species exhibits a decline in genetic diversity from East Asia to Europe [30], consistent with our results. This gradual reduction in genetic diversity along a cline could be explained by the “East Asian origin” of the gypsy moth proposed by many scholars, followed by the gradual spread through Central Asia to Europe and the artificial introduction of gypsy moths in North America in the early 19th century [27,28,30]. However, it has also been proposed that L. dispar originated in the early Pleistocene in the outer Caucasus, with a rapid early eastward spread (due to strong female flight ability), finally arriving in Japan, and later colonizing China via Russia and Mongolia [32]. We were unable to obtain reliable samples from Central Asia and therefore could not evaluate this hypothesis from the perspective of the nine microsatellite loci in this study. Based on AMOVA, we observed the majority of genetic variation in the 23 localities within populations, but we also detected genetic differentiation. Populations outside of China had relatively low genetic diversity, and differentiation may be explained by the separation of habitats and the low flight capacity of females. Our results further suggested that gene flow is higher in China (Table 3, Figure 3).

4.2. Other Subspecies Introduced to Eastern Eurasia

According to a review of potentially invasive Lymantria species by Pogue and Schaefer in 2007, L. dispar dispar is mainly distributed west of the Ural Mountains, in the Middle East, and in some islands in the Mediterranean, and is newly introduced in North America. L. dispar asiatica is mainly distributed east of the Ural Mountains to the Korean Peninsula and over two-thirds of China, among other places. L. dispar japonica is mainly found in Honshu Island, Shikoku Island, and parts of Kyushu Island, southwest of Hokkaido. However, it is unlikely that individuals with genotypic similarity to L. dispar dispar have emerged in Kuduer and Xifeng in China. Similar results have been obtained in a study of mitochondrial markers, that moths from Xifeng and Kuduer fall into the EGM clade in mitochondrial gene tree [34]. Samples collected in Kuduer were divided into a northern China cluster and a European cluster, and the Xifeng samples were classified into a southern China cluster and a European cluster (Figure 6). Different from the mitochondrial variation results of Zhao et al., the North American and European genotypes were closely related in the phylogenetic tree, while, based on nuclear genes, the gypsy moth genotypes in the Kuduer and Xifeng regions were more closely related to the gypsy moth genotypes in Europe than to those in North America. STRUCTURE implements different classification strategies for different loci and can be used to visualize the assignment of individual lineages to groups [57]. According to the results of genetic structure analysis, several individuals in Xifeng and Kuduer exhibited pure European cluster genotypes with the typical EGM, suggesting that it was recently introduced and has not yet fully interbred with local populations. The lack of private alleles in Europe and North America (Table 3) also suggests that the EGM spread in East Asia, mixing with local populations in other regions, leading to the widespread distribution of the EGM nuclear genotypes. A bottleneck test also indicated excessive heterozygosity caused by the recent introduction of EGMs in China (Table 3). The L. dispar dispar genotype was also found in Liaoning and Sichuan in China, further confirming the invasion of EGMs in China [35,75].

There is also evidence for a mixed genetic composition in Far East Asia. In an analysis of the population structure of gypsy moths in the boundary region, the peninsula, and adjacent areas of the Korean Peninsula, various L. dispar japonica genotypes in Japan were discovered in the southern part of the Korean Peninsula [31]. Wu et al. (2015) revealed a gradient of allelic dominance from L. dispar asiatica to L. dispar japonica from the Russian Far East to South Korea [30]. We obtained similar results in our analysis of microsatellite markers. Although the Primorsky Krai in Russia showed individuals assigned to the typical L. dispar asiatica group in northern China, there were genotypes of gene loci similar to L. dispar japonica. As an important port area in eastern Russia, we speculate that some individuals distinct from L. d. asiatica may have been introduced into the Russian Far East through the port.

4.3. Morphological Diversity

Biological diversity refers to both genetic and morphological diversity. Species in different environmental conditions can be differentiated at both the genetic and phenotypic levels, and species with low dispersal capacity might be more sensitive to local conditions [13]. Changes in traits under sexual selection in arthropods which are related to reproduction and affecting gene flow, have a major impact on species evolution [76], especially male genitalia in insects. Combined with analyses of insect morphological features and molecular genetic factors, comprehensive taxonomic analyses have provided key insights and enabled the discovery of species variation.

In a quantitative analysis of male genitalia of gypsy moths from different regions, we did not find clear distinctions among subspecies overall, other than the single subspecies in Honshu Island, which was significantly different from other subspecies with respect to genitalia. Genitalia of gypsy moths in the Russian Far East showed similarities to those of the Eurasia gypsy moth and the Japanese unique subspecies. This also indirectly illustrates the mixed composition of the gypsy moth population in the Russian Far East. Yunnan and Zunyi in the south showed separation from gypsy moth populations at higher latitudes in China with respect to genitalia. Additionally, male genitalia of gypsy moths in the United States, believed to have been imported from Europe, also showed differences from those of gypsy moths in Europe. Apparently, the weak flying ability and poor mobility of larvae make L. dispar dispar highly sensitive to differences in environmental conditions, resulting in gradual differentiation from the gypsy moth in Europe after its introduction to North America [77,78,79].

4.4. Unique Clusters in Southern China

Gypsy moths in China are assumed to belong to a single lineage of L. dispar asiatica, and most studies are limited to samples from northern China. However, we disproved this assumption based on analyses of samples from across China. In addition to individuals with EGM genotypes in China, we detected differentation between gypsy moths in southern and northern China. Despite similar levels of genetic diversity, Chengdu and Yunnan had the most private alleles, and the F_ST value for comparison between south of Anhui and north of Anhui exceeded 0.05, indicating moderate or above-average genetic differentiation. This was not observed between populations in the north, and even gypsy moths in Yunnan might have a high degree of genetic differentiation. The results of a STRUCTURE analysis (K = 5) clearly showed the difference between gypsy moths in the southern part of Anhui and those in the northern part of China (Figure 6). Anhui is at the junction of northern and southern China, and the genotypes of gypsy moths in this area illustrate a degree of north–south lineage mixing.

Little is known about the current population structure of gypsy moths in southern China. Zhao et al. (2019) revealed from mtDNA analyses that sister groups of typical L. dispar asiatica in northern China are present in southern China. We confirmed this result using nuclear genes. Based on haplotype differences, differentiation in the Zunyi and Dayi samples may be explained by the warm climate in the south, as the two regions serve as a glacial refuge for gypsy moths [33]. The strong geographical isolation between the shelters and secondary exposure to post-glacial conditions after the temperature rise might explain the genetic differentiation of gypsy moths in the south, but further studies are needed to evaluate this hypothesis. In addition, we analyzed the north–south geographical boundary and found that the fuzzy boundary between groups may correspond with the Qinling–Huaihe line of China. The Qinling Mountain–Huaihe River line is the traditional climate divide between temperate and sub-arid North China and subtropical and humid South China [80]. It includes mountains and rivers as geographical barriers. Although gypsy moths in China have stronger flight ability than EGMs, the potential flight distance is likely to be insufficient to traverse the geographic barrier, preventing migration between the two temperature zones [81,82]. Intraspecific variation in subtropical and temperate regions, and differentiation between the two temperature zones, are well-established [83,84,85]. In an analysis of male genitalia, we found that individuals in Zunyi and Yunnan were clearly morphologically distinct from individuals in other populations. We supposed that in addition to the role of the interglacial period, gypsy moths have inhabited different climatic, geographical, and host environments for a long time in south and north in China. The weak migration ability of the species may make it sensitive to local environmental conditions, resulting in distinct patterns of genetic differentiation.

5. Conclusions

In summary, we found evidence for both EGMs and Japanese gypsy moths in East Asia, suggesting the potential for extensive population mixing. In East Asia, we detected high genetic diversity, indicating that the gypsy moths in this region are more adapted to the environment through the variation than the gypsy moths in Europe and North America. Their strong diffusion ability and wide host range might enable migration to non-native areas, resulting in ecological damage and economic losses. We did not detect AGMs in North America, indicating that strict quarantine measures successfully prevented the establishment of AGMs via goods and vessels. A group of Central and South China regions with moderate or high genetic differentiation from typical L. dispar asiatica showed a high degree of genetic and morphological variation over a wide range of latitudes and longitudes. In future studies, additional regions, especially those in port areas, should be evaluated using more extensive methodologies to explain the spread of the EGM to China and to clarify whether the genetic differentiation in southern China has reached the level of subspecies differentiation.