Mitochondrial DNA markers reveal high genetic diversity and strong genetic differentiation in populations of Dendrolimus kikuchii Matsumura (Lepidoptera: Lasiocampidae)

Dendrolimus kikuchii Matsumura, 1927 is a serious forest pest causing great damage to coniferous trees in China. Despite its economic importance, the population genetics of this pest are poorly known. We used three mitochondrial genes (COI, COII and Cytb) to investigate the genetic diversity and genetic differentiation of 15 populations collected from the main distribution regions of D. kikuchii in China. Populations show high haplotype and nucleotide diversity. Haplotype network and phylogenetic analysis divides the populations into three major clades, the central and southeastern China (CC+SEC) clade, the eastern China (EC) clade, and the southwestern China (SWC) clade. Populations collected from adjacent localities share the same clade, which is consistent with the strong relationship of isolation by distance (r = 0.74824, P = 0.00001). AMOVA analysis indicated that the major portion of this molecular genetic variation is found among the three groups of CC+SEC, EC and SWC (61.26%). Of 105 pairwise FST comparisons, 93 show high genetic differentiation. Populations of Puer (PE), Yangshuo (YS) and Leishan (LS) are separated from other populations by a larger genetic distance. Distributions of pairwise differences obtained with single and combined gene data from the overall populations are multimodal, suggesting these populations had no prior population expansion in southern China. The nonsignificant neutral test on the basis of Tajima’ D and Fu’s Fs, and the lack of a star-shaped haplotype network together with the multiple haplotypes support this hypothesis. Pleistocene climatic fluctuations, combined with the host specificity to Pinus species, made these regions of south China into a refuge for D. kikuchii. The high level of population genetic structuring is related to their weak flight capacity, their variations of life history and the geographic distance among populations.


Introduction
Simao pine moth, Dendrolimus kikuchii Matsumura, 1927, is a destructive forest pest with an extensive range across southern China. Larvae attack various coniferous trees and regular outbreaks occur. These outbreaks are mainly on Pinus langhianensis Chev, P. yunnanensis Franch, P. massoniana Lambert, P. armandi Franch and Keteleeria evelyniana Mast. [1,2]. D. kikuchii larvae consume, on average, 7486.6 mm of P. langhianensis pine needles to complete their development [3]. During outbreaks, a large amount of pine needles are consumed, giving the appearance of forest fire damage. Larval damage from a large infestation of D. kikuchii can reduce the yield of timber, resin and cones, and affect the growth rate of pines, even resulting in tree death [1,2]. The Simao pine moth shows differences in ecological preferences among areas in China. One generation of D. kikuchii occurs in areas with short periods of optimal environment, such as in the middle region of Yunnan and Guizhou provinces. In Zhejiang, Fujian, and the southwest region of Yunnan provinces, two generations occur. In different counties of Yunnan province, such as Jingdong and Anning, the life histories of D. kikuchii also show significant variation [1].
Genetic diversity and genetic structure in insects can be affected by many factors, such as host plant species, chemical controls, geographic distance and geographic barriers [4][5][6][7][8][9][10][11]. For some lepidopterous species, the genetic diversity and genetic structure are related to their migration capacity and number of generations [12][13][14][15]. Despite the economic and landscape threats of the D. kikuchii to pine trees in the southern parts of China and the need to establish control strategies for this pest, it has been unclear whether analogous associations between these factors and genotype patterns may be present among D. kikuchii populations.
Pleistocene climatic fluctuations are thought to have a great effect on this species' distribution and historical demography [16]. South China was considered as a key area of some refugia for many relict rare species during the Pleistocene glaciation [17,18]. Zhang et al [19] studied the geohistory of Dendrolimus punctatus, a sympatric species of D. kikuchii, based on these geological events. They found that D. punctatus settled down in south China with the spread of masson pine during the Pleistocene. However, to our knowledge, there is no published report on the population history of D. kikuchii, and also no studies on the population history of Dendrolimus species based on molecular data.
In this study, we use three mitochondrial genes to (i) investigate the genetic diversity and genetic differentiation of 15 populations collected from the main distribution regions of D. kikuchii in southern China, and infer the demographic history of this pest, and (ii) to test the hypothesis that geographical isolation and biological characters, such as life history, are significant factors underlying genetic variation in Chinese D. kikuchii.

Ethics statement
There is no endangered or protected species involved in this study. No specific permissions were required for the described field studies for this widespread forest pest. We confirm that the locations are not privately owned or otherwise protected.

Sampling
A total of 182 individuals were collected from 15 locations during 2013 to 2016 within the D. kikuchii China distribution range (Table 1 and Fig 1). Different instar larvae and adult moths were sampled. For larvae, only one individual was collected per tree. For adult moths, pheromone traps were used in the pine forest (>1 ha) and only one moth from each trap was sampled in order to avoid sampling errors. All sampled insects were immersed in absolute ethyl alcohol, and then stored at -20˚C prior to genetic analysis.

DNA extraction and amplification
Genomic DNA was extracted from the last proleg of the larvae or from any leg of the adult using the DNeasy Tissue Kit (QIAGEN, Hilden, Germany). Extraction was performed according to the bench protocol for animal tissues. Mitochondrial COI (LCO1490, HCO2198, [20]), COII (Eva, Patrick, [21]), and Cytb (generated for codling moth, unnamed [22]) were selected for use in this study. Polymerase chain reactions (PCR) were performed using a S1000 Thermal Cycler (BIO-RAD, Hercules, CA, USA) in a total volume of 20 μl, containing 10 μl 2×PCR Super Master Mix (Biotool, Shanghai, China), 0.25 μM of each primer, and 1 μl genomic DNA (10-30 ng/μl). PCR amplification was employed with denaturation at 95˚C for 10 min, followed by 40 amplification cycles consisting of 95˚C for 30 s, primer-specific annealing temperature of 53˚C (COI), 52˚C (COII and Cytb) for 1 min, 72˚C for 45 s, and then a final step at 72˚C for 10 min. Amplified products were purified and sequenced by Tianyi Huiyuan Biotechnology Co., Ltd.

Data analysis
The sequences were preliminarily aligned using the CLUSTAL X program [23]. Sequences of COI (646), COII (675), Cytb (700) of D. kikuchii were deposited in the NCBI GenBank (Gen-Bank accession numbers: COI, MF155667-MF155697; COII, MF155698-MF155737; Cytb, MF155738-MF155760) (data in S1-S3 Text). Multiple sequences of COI, COII and Cytb were concatenated to yield a total length of 2021 bp. The haplotype network of D. kikuchii was analyzed using a median-joining algorithm in the program Network 4.6 [24]. A neighbour-joining (NJ) tree was built using NJ tree subroutine in NEIGHBOUR within PHYLIP 3.5 [25], and the parameters were expanded 1,000 times. The CONSENSE subroutine within PHYLIP was then applied to generate a consensus NJ tree that provided bootstrap support at each node. The tree was visualized using treeview version 1.6.6 software. DnaSP 5.0 [26] were performed to calculate number of polymorphic sites (S), number of haplotypes (H), haplotype diversity (Hd), nucleotide diversity (Pi), Tajima's D (D), and Fu's Fs (Fs). Analysis of molecular variance (AMOVA) was performed using the ARLEQUIN version 3.5 based on the combinations of the three gene sequences [27], along with calculating pair fixation indices (FST). The pairwise genetic distances were calculated by MEGA 6.0 [28] based on the Kimura-2-parameter model [29]. Referring to the criterion for genetic differentiation by Wright (1978) [30], we defined genetic differentiation as low for FST<0.05, moderate for 0.05<FST<0.15, high for 0.15<FST<0.25, and very high for FST>0.25 [31]. In order to test isolation by distance (IBD), the matrices of genetic distance FST/(1-FST) and the geographic distance (ln) between all 15 sampling populations were compared using the Mantel test with 10,000 permutations [32]. This analysis was performed using the ZT software package [33]. For examining demographic history, the distribution of pairwise differences between individual sequences was analyzed by means of mismatch distribution analysis using DnaSP 5.0 [26]. The formula, Tau = 2 ut was used to detect the time of population size changes [34]. The nucleotide substitution rate in mitochondrial DNA was 2.3% per million years (MY) as suggested in Knowles et al. (2000) [35].  Table 1; the gray region represents the geographic distribution of five D. kikuchii populations in Anhui. https://doi.org/10.1371/journal.pone.0179706.g001 Population genetics of Dendrolimus kikuchii in China

Genetic diversity
For concatenated sequences, the haplotype diversity ranged from 0 to 1 with a average of 0.940, while the nucleotide diversity ranged from 0 to 0.00251 with a average of 0.01451 ( Table 2). All populations displayed large numbers of mitochondrial haplotypes, with a total 31 haplotypes obtained for COI, 40 haplotypes for COII, and 23 haplotypes for Cytb. Seven, six and four common haplotypes were shared for COI, Cytb and COII respectively (Fig 2). The median-joining network demonstrated a high genetic diversity for populations of D. kikuchii in China. Of 73 examined haplotypes, 68 were unique. Most populations lacked a common haplotype.

Genetic differentiation
The Median-Joining network of the haplotypes can be divided into three major clades (Clades CC+SEC, EC and SWC) (Fig 3). All haplotypes in Clade CC+SEC were obtained from samples from central China and southeastern China including Hunan, Jiangxi, and Fujian provinces, while all haplotypes in Clade EC were from eastern China including Anhui and Zhejiang provinces. The haplotypes in Clade SWC included samples from southwestern China and covered Yunnan, Guizhou and Guangxi provinces. Overall, populations collected from adjacent localities or the same province shared the same clade. However, the Median-Joining network of the haplotypes based on single gene of COI and Cytb was not in line with the result from the combined genes. This may due to less nucleotide variation in both single genes. Phylogenetic reconstruction using the Kimura-2-parameter resulted in a consensus NJ tree with comparatively higher bootstrap values. The populations were divided into three major clusters, CC+SEC, EC and SWC, which comply well with the results from the haplotype network (Fig 4).   (Table 3). This is consistent with the results of the clustering analysis based on a Median-Joining network.  (Table 3).
AMOVA results indicate that the major portion of the molecular genetic variation is found among groups (61.26%). Exact tests showed a significant genetic variance on all three levels (P<0.001) ( Table 4).
The Mantel test for the 15 populations revealed a positive correlation between genetic distances and geographic distances (r = 0.74824, P = 0.00001), suggesting that isolation by distance had a limiting effect on gene flow.

Mismatch distribution
The results of the combined gene analysis show that Tajima' s D values are significantly positive with a value of 0.24250, but are not significant in most specific populations (P<0.05 in HS, PE and HY, P>0.05 in the rest of the populations). Fu's F statistic was significantly negative with a value of -5.132 (P>0.1) ( Table 5). At the population level, HS, HY and PE populations have negative and significant Tajima' s D and Fu's Fs values, whereas HS, HY and PE populations show a bimodal distributions, suggesting that expansion events were not detected with this analysis (data not shown). Distributions of pairwise differences (mismatch distributions) obtained with the single and combined gene data from the overall populations were multimodal, suggesting that the populations of D. kikuchii in southern China did not experience population expansion (Fig 5). A nonsignificant neutral test based on Tajima' D and Fu's Fs support this interpretation. The time of reaching a stable population size are estimated to be 26500 (COI), 21600 (COII), and 11800 (Cytb) years ago.

Discussion
Using three mitochondrial genes, we investigated the genetic diversity and structure of 182 individuals of 15 D. kikuchii populations sampled throughout their main areas of distribution in China. The results show a high genetic diversity and high level of genetic structuring of D. kikuchii in the sampled areas.
All populations displayed large numbers of mitochondrial haplotypes, with a total 31 haplotypes obtained for COI, 40 haplotypes for COII, and 23 haplotypes for Cytb, of which seven, six and four were common haplotypes shared respectively. Of 73 examined haplotypes based on combined genes, 68 were unique and did not share the same ancestral haplotype. High numbers of private haplotypes and lack of ancestral haplotype suggest that D. kikuchii could be a species native to China. In a previous study, Zheng et al. [36] indicated that these two features can be considered as indicators for native species, especially for native species with a low   dispersal capacity, such as Grapholita molesta Busck (Lepidoptera: Tortricidae) and Chilo suppressalis (Walker) (Lepidoptera: Pyralidae). Although the origin of D. kikuchii remains uncertain, the higher genetic diversity we observed in the present study supports our assumption that China is part of the original range of D. kikuchii. Based on the results of three mitochondrial genes, including AMOVA analysis, haplotype network and phylogenetic analysis, we conclude that D. kikuchii populations were geographically structured in three regions: eastern China, southwestern China, and central China along with southeastern China. The sampling localities of southwestern China are isolated by the Wumeng, Leigong, and Fanjingshan Mountains, while sampling localities of central China are isolated by the Hengshan and Dabieshan Mountains, which have acted as substantial barriers to gene flow. Populations from central China and southeastern China shared common haplotypes; this is in accordance with the results of phylogenetic analysis. Although there is a geographical barrier formed by Wuyishan Mountain between these two regions, a transportation network increases the gene flow among D. kikuchii populations. In Lepidoptera species, dispersal patterns also influence genetic variation [37,38]. D. kikuchii is generally regarded as a sedentary species based upon previous studies of its flight capacity, with a sphere of activity extending over 10~20 km [1,3]. The weak flight capacity of D. kikuchii can reduce gene flow among populations. The IBD relationship (r = 0.74824, P = 0.00001) in the present study supports this hypothesis. Similar trends have been identified for many sedentary species, such as Chilo suppressalis (Walker) [15] and Carposina sasakii Matsumura [18]. The same result was reported by Weng et al. [39]. They found that the genetic variation of D. punctatus populations that ranged over five adjacent regions in Zhejiang province were low using the ISSR-PCR marker.
Analysis of genetic distance indicate that populations of Puer (PE), Yangshuo (YS) and Leishan (LS) are separated from other populations. In fact, D. kikuchii caterpillars in Puer and Yangshuo populations turn into adults one moth earlier than those of other populations. And in Guizhou, only one generation occurs, while two generations occur in the rest of the provinces [1]. These variations in life history may contribute to the significant genetic differentiation. Moreover, the areas of distribution of D. kikuchii in China occur across very complex topography (tall mountains, plain and basin), different climates (temperate and tropical climates), different agricultural landscapes and forest types. A high level of population differentiation presents a high potential for adaptation to different environmental conditions [40] and high reproductive rates [41], which allow the moth to form locally differentiated populations. Similar results have been found for Chilo suppressalis [15] and Reticulitermes chinensis [42].
Distributions of pairwise differences (mismatch distributions) obtained with COI, COII, Cytb, and combined gene data from the overall populations are multimodal, suggesting that the populations of D. kikuchii did not experience population expansion. The lack of a star shape for the haplotype network together with the existence of multiple haplotypes support such a hypothesis. This also explains why there are so many private haplotypes in different geographical populations. The times when population histories stabilized are estimated to be 26500 (COI), 21600 (COII), and 11800 (Cytb) year ago. This range is within the late Pleistocene, which was characterized by climatic oscillations between warm and cold periods [43]. Over 42 thousand years ago, Pinus massoniana, the favorite host plant for D. kikuchii, gradually spread into the warmer south China, a condition they would prefer [19]. These events impacted the distribution of D. kikuchii, with these regions of south China gradually becoming the refuge for D. kikuchii. During the late Pleistocene, some regions of south China were not covered by large ice sheets [44]. In addition, south China comprises a mountainous mosaic area and has the potential to host microclimatic zones that are probably capable of providing some key refuges for many relict species [45]. Considering the lack of an ancestral haplotype and strong isolation-by-distance relationships of this species, we can conclude that D. kikuchii in south China has arisen in separate refuges and experienced parallel evolutions. Similar results have been found for grasshoppers [46], and the blue manakin [47].
Future population genetic research on D. kikuchii in China should cover a larger area and larger number of sampled individuals, as well as use nuclear genes as markers to provide a greater understanding of genetic structure. Additional research needs to be done to detail the geographic origins of D. kikuchii in China and its spread through China.

Conclusion
Using three mitochondrial genes, we investigated the genetic diversity and structure of 182 individuals of D. kikuchii sampled throughout its main distribution areas in China. The results show high genetic diversity and a high level of genetic structuring of D. kikuchii in sampled areas. The high level of population genetic structuring is related to the weak flight capacity of the D. kikuchii, variations in its life history and the geographic distance among populations. Distributions of pairwise differences (mismatch distributions) obtained with COI, COII, Cytb, and combined gene data indicate the populations of D. kikuchii in southern China did not experience population expansion. These genetic data not only provide us with an understanding of population genetics for such a secondary species, but also provide guidance for pest management strategies.