To elucidate the genetic population structure of two coastal weevils, Scepticus griseus and S. tigrinus, we conducted molecular phylogenetic analyses of the mitochondrial DNA cytochrome c oxidase subunit I (COI) region (1308 bp) and cytochrome c oxidase subunit II (COII) region (584 bp). A total of 650 individuals (S. griseus, 444 individuals; S. tigrinus, 206 individuals) were obtained from 64 sites. The haplotype networks of both species showed three major lineages with roughly regional distribution. However, the two species show quite different genetic structures; S. griseus has a complicated structure while that of S. tigrinus is simple. We hypothesize that the genetic structure of each of these two weevil species reflects climatic oscillations during the Pleistocene, and the differences in genetic structure between S. griseus and S. tigrinus may represent a unique evolutionary history scenario in each species.
INTRODUCTION
The distribution of animals along seashores can be treated as one-dimensional (Pielou, 1979) and, consequently, analyses of population structure and relationships between populations are more easily conducted with results that are more simply interpretable compared with three-dimensional distributions. Numerous studies of population structure in insects and isopods adapted to seashore environments have been conducted (Satoh et al., 2004; Hojito et al., 2010; Kudo et al., 2012; Niikura et al., 2015), and while their distributions show similarities (e.g., the divergence between the Pacific Ocean and the Sea of Japan sides), these studies also show that each species has a unique population structure. For example, Niikura et al. (2015) showed that the genetic structure of Tylos granuliferus on the main islands of Japan is strongly related to oceanic currents. In contrast, several coastal species showed population structure with geohistorical influence (Satoh et al., 2004; Itani, 2000). Therefore, the population structure of coastal species is likely affected by either oceanic currents or geohistorical factors, which show noncurrent-related distribution patterns.
Apterous curculionid weevils, Scepticus griseus (Roelofs, 1873) (= S. uniformis Kôno, 1930) and S. tigrinus (Roelofs, 1873), belong to tribe Tanymecini of subfamily Entiminae. Both inhabit sandy beaches, but with rough separation in Japan with S. griseus is distributed on the Pacific Ocean side and S. tigrinus on the Sea of Japan side of the archipelago (Hokkaido, Honshu, Shikoku, and Kyushu) (Morimoto, 1962; Sawada, 2008; Morimoto et al., 2015). These two species are very similar to each other in morphological appearance, but molecular analyses using the mitochondrial cytochrome c oxidase subunit I (COI) region clearly show that the two are separate species (Yamashita et al., 2015). As both species have limited migration ability due to the flightlessness of the species and genetic differences were recognized at the population level in both species (Yamashita et al., 2015), they are expected to have complex genetic structures. Interestingly, the two species show color variation irrespective of geographical trends or localities, and Sawada (2008) suggested that gene flow occurs between populations of the two species by current activity.
Therefore, we can expect that the distributions of these two weevils show current- or noncurrent-related distributional patterns. In order to assess the population structure and phylogeny of the two weevils and elucidate the historical events responsible for their present distributions, we analyzed COI and cytochrome c oxidase subunit II (COII) gene sequences from 64 localities.
Table 1.
Sampling sites, number of specimens (N), haplotype information, haplotype diversity (Hd), nucleotide diversity (Pi), Tajima's D and Fu's Fs of S. griseus and S. tigrinus. Asterisk (*) indicates the unique haplotype excluded from genetic analyses excepting haplotype network analyses. Double asterisk (**) indicates p < 0.05.
Continued
Continued
MATERIALS AND METHODS
We collected specimens of S. griseus and S. tigrinus from September 2010 to June 2015 at 64 localities in Japan and Korea (Table 1 and Fig. 1). A total of 650 individuals (S. griseus, 444 individuals; S. tigrinus, 206 individuals) were sampled from under or on weeds along the seashore. All specimens were stored in 99.5% ethanol until DNA extraction.
To prepare specimens for molecular analyses, all legs on the right side of the body were removed, and DNA was extracted using a DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany). The legs were homogenized in 180 µl ATL buffer with 20 µl proteinase K and incubated at 55°C for 1 h. After incubation, total genomic DNA was extracted following manufacturer instructions.
To characterize the population structures of the two species, we analyzed the COI region (S. griseus, 508 bp; S. tigrinus, 503 bp) in 444 individuals of S. griseus and 206 individuals of S. tigrinus. For phylogenetic analysis within populations, additional regions were sequenced for 23 S. griseus and 19 S. tigrinus of the individuals randomly selected from among the initial sample: another COI region (S. griseus, 800 bp; S. tigrinus, 805 bp) (Fig. 2) and COII region (622 bp).
We amplified the COI and COII genes by polymerase chain reaction (PCR) using the primers shown in Table 2. Because the target gene region of S. griseus did not amplify well, we designed a new forward primer, COI-F, which bound to a site inside the primer region that described by Folmer et al. (1994) (Table 2). PCR thermal cycling conditions for the COI-1 set (LCO1490 and HCO2198) and COI-2 set (COI-F and HCO2198) were 94.0°C for 2 min, 35 cycles at 94.0°C for 30 sec, 49.0°C for 1 min, and 72.0°C for 1 min, with a final extension at 72.0°C for 5 min. PCR thermal cycling conditions for the COI-3 set (C1-J-1718 and TL2-N-3014) and COI-4 set (C1-J-2183 and TL2-N-3014) were 94.0°C for 5 min, 40 cycles at 94.0°C for 30 sec, 47.0°C for 30 sec, and 72.0°C for 1 min, with a final extension at 72.0°C for 5 min. To amplify the COII region, we used primers TL2-J-3037 and C2-N-3661 (Table 2) with PCR thermal cycling conditions of 94.0°C for 5 min, 35 cycles at 94.0°C for 30 sec, 49.0°C for 30 sec, and 72.0°C for 30 sec, with a final extension at 72.0°C for 5 min. PCR was performed in a reaction volume as was described by Yamashita et al. (2015). The amplified DNA was sequenced directly by an automated method using the BigDye® Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, USA) on an Applied Biosystems® 3730xl DNA Analyzer. These sequence data have been deposited in GenBank (accession numbers: LC108859–LC108865, LC108867–LC108871, LC108879–LC108885, LC108887–LC108889, LC108892, LC108894, LC108896–LC108900, LC108902–LC108903, LC108906, LC108908, LC108911, LC108913, LC108915, LC108917–LC108918, LC108920, LC108923, LC108925, LC108931, LC108960, LC108969, LC108971–LC108979, LC108982– LC108984, LC108986–LC108988, and LC109324–LC109620).
Fig. 1.
Sampling localities for S. griseus and S. tigrinus, with population numbers (referenced in Table 1). Closed circles indicate populations included in phylogenetic analysis.
Table 2.
Primers used for the amplification of the fragments of the mitochondrial DNA COI and COII regions.
Sequences were aligned using Muscle (Edgar, 2004a, b) in MEGA version 5 software (Tamura et al., 2011). To determine whether the congruence between the mitochondrial DNA COI and COII regions was significant, the incongruence length difference (ILD) test (Farris et al., 1995) was conducted by excluding all invariant sites and using the partition homogeneity test in PAUP ver. 4.0 beta10 (Swofford, 2002) with 1000 replicates. As there were no significant differences (P > 0.05), we also used these two genes (COI and COII) as a combined dataset (1892 bp) for molecular phylogenetic analyses.
Fig. 3.
Median-joining (MJ) networks of COI haplotypes detected in S. griseus. Circle size indicates sample size. Colors of circle reflect each locality. Missing haplotypes are indicated by black circles.
Phylogenetic analyses were performed by the maximum likelihood (ML; Felsenstein, 1981) method, using MEGA software. ML analysis was conducted using the general time reversible model (GTR) based on the log likelihood criteria. To estimate the confidence probability for each interior branch, the bootstrap method (Felsenstein, 1985) was performed with 1000 replications. In the phylogenetic analysis, Scepticus insularis Roelofs, 1873 (= S. griseus auct., nec Roelofs, 1873) was used as the outgroup. All gaps were treated as missing data with the complete-deletion option. We then constructed a ML tree on the basis of the analysis. Estimates of the divergence times of S. griseus and S. tigrinus were calculated with MEGA version 5. We used a rate of 2.3%/Ma applied a “standard” mitochondrial DNA clock for insects (Brower, 1994).
DnaSP v5 (Librado and Rozas, 2009) was used to search for haplotype variation, and NETWORK 4.5.1.0 ( http://www.fluxus-engineering.com) was used to construct a median-joining (MJ) network (Bandelt et al., 1999) to describe relationships of the haplotypes.
We conducted population genetic analyses using sequences from the 40 populations for which more than five individuals were available. In order to estimate the extent of genetic diversity in each population, we calculated nucleotide diversity (π) (Nei, 1987) with DnaSP v5. Tajima's D (Tajima, 1989a, b) and Fu's Fs statistics (Fu, 1997) were calculated to infer the extent of population expansion with ARLEQUIN 3.5. To estimate the extent of divergence between populations, we calculated the net nucleotide substitutions between populations (d A) (Nei, 1987) with MEGA version 5. Finally, to estimate the relative degree of genetic differentiation, we calculated pairwise FST values among populations with more than 10 individuals using ARLEQUIN 3.5, and we tested the significance of the FST values following 1000 permutations in ARLEQUIN 3.5. We performed an analysis of molecular variance (AMOVA; Excoffier et al., 1992) in ARLEQUIN 3.5 (Excoffier and Lischer, 2010) with significance tests conducted using a non-parametric permutation approach (1023 permutations performed) to examine the genetic structure of populations.
Voucher specimens from this study are preserved in the Laboratory of Entomology, Tokyo University of Agriculture (TUA), Atsugi, Japan.
RESULTS
Haplotypes and haplotype relationships
Analysis of partial nucleotide sequences of the mitochondrial COI coding region (S. griseus: 508 bp; S. tigrinus: 503 bp) showed a total of 214 haplotypes among 444 Individuals collected from 47 sites in S. griseus, and a total of 83 haplotypes among 206 individuals collected from 22 sites in S. tigrinus, respectively (Table 1 and Fig. 1). Among the 214 haplotypes, 199 were unique to single sampling sites and the others were shared by two to six populations in S. griseus, and among the 83 haplotypes, 75 were unique to single sampling sites and the others were shared by 10 populations in S. tigrinus.
Figures 3 and 4 show haplotype networks of each species constructed by the MJ method. Three lineages are roughly recognized in both S. griseus and S. tigrinus, and are designated here as lineage A, lineage B and lineage C in S. griseus, and lineage D, lineage E and lineage F in S. tigrinus. None of the lineages of S. griseus showed star topologies. The dA values between lineages A, B and C ranged from 0.042 to 0.073 (Table 3). On the other hands, in S. tigrinus, lineage D showed a star-like topology with one haplotype (Hap_T5) shared by 10 populations (Table 1) at the center, excluding one exception (Hap_T10). Lineages E and F each displayed a scattered network, and lineage E was divided into three sub-lineages with at least seven nucleotide substitutions between the sub-lineages. The dA values between lineages D, E and F ranged from 0.029 to 0.043 (Table 4).
The geographic distributions of the three lineages of S. griseus were nearly mutually exclusive (Fig. 5). Haplotypes of lineage A occurred in Kanto, which includes Fukushima Prefecture, the Izu Peninsula and the Izu Islands (sites 1–16 in Table 1 and Fig. 1); haplotypes of lineage B were found in Tokai, Kii and Setouchi (sites 17–26 in Table 1 and Fig. 1); and haplotypes of lineage C were found in the Ryukyus, Kyushu, westernmost Honshu (Yamaguchi and Shimane Prefectures), southern Shikoku (Kochi Prefecture) and Korea (sites 27–47 in Table 1 and Fig. 1). The populations were thus clearly distinguishable by their haplotypes. Only localities 11 and 16 had haplotypes that were mixes of lineages A and B. The geographic distributions of the three lineages of S. tigrinus were also mutually exclusive, but strictly (Fig. 5). Haplotypes of lineage D occurred in Hokkaido and Tohoku (including Sado Island and Niigata Prefecture; sites 48–56 in Table 1 and Fig. 1); haplotypes of lineage E were found in Ishikawa, Fukui and Kyoto Prefectures (sites 57–59 in Table 1 and Fig. 1); and haplotypes of lineage F were found in San'in and northern Kyushu (sites 29, 31–34 and 60–64 in Table 1 and Fig. 1). The populations were thus clearly distinguishable by their haplotypes. Only localities 55 and 58 had haplotypes that were mixes of lineages D and E.
Fig. 4.
Median-joining (MJ) networks for the COI haplotypes detected in S. tigrinus. Circle size indicates sample size. Colors of circle reflect each locality. Missing haplotypes are indicated by black circles.
Table 3.
Pairwise FST (below diagonal) and dA (above diagonal) values for lineages of S. griseus.
Table 4.
Pairwise FST (below diagonal) and dA (above diagonal) values for lineages of S. tigrinus.
Table 5.
Haplotype diversity (Hd), nucleotide diversity (Pi), Tajima's D, and Fu's Fs for six lineages of S. griseus and S. tigrinus. Asterisk (*) indicates p < 0.05.
Genetic diversity and differentiation
Haplotype diversity of S. griseus was high both within species and within lineages (Tables 1 and 5). Nucleotide diversity within lineage C was higher than in lineages A and B. Pairwise FST values between lineages and between populations in each lineage (Tables 3 and 6) were significant (P < 0.05). AM OVA identified the largest genetic variance as being among lineages followed by the variance among populations (Table 7). Fu's Fs showed significant negative values in only eight populations (sites 6, 9, 20, 21, 22, 35, 42 and 43 in Table 1). Tajima's D showed significant negative values in only one population (site 9 in Table 1).
In S. tigrinus, the haplotype and nucleotide diversities of lineage D were lower within lineages than those of lineages E and F (Table 5). Pairwise FST values between lineages and between populations in each lineage (Tables 4 and 8) were significant. The values between populations in lineage D were lower than those in lineages E and F. AMOVA identified the largest genetic variance as being among lineages followed by the variance among populations (Table 7).
Fu's Fs showed significant negative values in lineage D, in lineage F, in four populations within lineage D, and in two populations within lineage F (sites 48, 50, 53, 54, 60 and 61 in Table 1). Tajima's D showed significant negative values in lineage D, in lineage F, in four populations within lineage D, and in one population within lineage F (sites 50, 52, 53, 54 and 60 in Table 1).
Estimation of divergence times
The phylogenetic tree and divergence times of lineages of each species are shown in Fig. 6. We described the divergence of lineages A + B and lineage C as node 1, lineage A and lineage B as node 2, lineages D + E and lineage F as node 3, and lineage D and lineage E as node 4. Node 1 was estimated to occur at 1.81 Ma, node 2 was at 0.63 Ma, node 3 was at 0.93 Ma, and node 4 was at 0.45 Ma. The estimated divergence times of the six lineages (1.81–0.45 Ma) was during the Early and Middle Pleistocene. According to evidence of climatic oscillation includes the Marine Isotope Stage (MIS), which serves as a barometer for sea level (Lisiecki and Raymo, 2005), 1.81 Ma and 0.63 Ma corresponded with MIS 64 and 15 in S. griseus, and 0.93 Ma and 0.45 Ma corresponded with MIS 24 and 12 in S. tigrinus (uneven numbers of MIS: interglacial periods; even numbers of MIS: glacial periods). The divergence times for S. tigrinus are later than for S. griseus.
Table 6.
Pairwise FST (below diagonal) and dA (above diagonal) values among 24 populations of S. griseus. Asterisk (*) indicates insignificant FST values (p > 0.05).
Table 7.
Analysis of molecular variance (AMOVA) for genetic structure within and among populations.
Table 8.
Pairwise FST (below diagonal) and dA (above diagonal) values among 15 populations of S. tigrinus.
DISCUSSION
Biogeographic pattern of Scepticus griseus
The results from the haplotype network show the following regarding the lineages [lineage A, lineage B (excluding Mie populations) and lineage C (Kyushu only)]: (1) a few haplotypes are shared between populations; (2) there are no major haplotype; (3) a number of haplotypes are connected by one step; (4) there are few missing haplotypes; and (5) unique haplotype networks are not found in populations. Thus, it can be inferred that the extinction of populations was rare and population size was relatively stable because Fu's Fs and Tajima's D showed significantly negative values in only eight (Fu's Fs) and one (Tajima's D) populations out of 40 populations (Table 1). Furthermore, the presence of various minor lineages in many populations in the haplotype network suggests that isolation of the respective populations and gene flow between populations was repeated over a short period. In addition, high haplotype and nucleotide diversities were recognized. Taken together, genetic diversities in the respective lineages have been maintained in this species.
The phylogenetic tree in Fig. 6 shows that genetic differentiation of node 1 (lineage A + B and lineage C) and node 2 (lineage A and lineage B) were estimated to occur during the Early and Middle Pleistocene when repeated climatic oscillation occurred. According to MIS proposed by Lisiecki and Raymo (2005), the genetic differentiation of node 1 occurred during the glacial period when climate was cold. While, that of node 2 was occurred during interglacial period when a marine transgression occurred. Repeated climatic oscillation caused multiple regressions and transgressions of the coast line in repeated glacial and interglacial periods. As S. griseus lives along seashores, its habitat was probably influenced by the regressions and transgressions; that is, population size of the species was drastically changed. In all, the lineages and populations were isolated, and gene flow was unlikely to have occurred between the populations during the interglacial periods. In contrast, the majority of the populations were in contact with each other through gene flow in the glacial periods.
The relatively high FST and dA values indicate the presence of great divergences between populations in each lineage, and it is likely that gene flow between the populations has been limited since the last glacial period due to habitat segregation. This is in contrast to the interglacial period when populations became isolated that was interspersed with gene flow between populations during glacial periods. The flightless weevil is presumed to always have had access to appropriate habitat, such as continuous sandy shores. Furthermore, extinction of populations seldom occurs due to the relatively warm climate influenced by a warm current (Kuroshio current) along the Pacific Ocean side. For these reasons, the genetic structure of this species is considered to have become complicated. High genetic divergences between the lineages is assumed to have developed through the prevention of gene flow by geological gaps such as large rivers and extended stretches of rocky shores.
Fig. 6.
(A) Phylogenetic relationships of S. griseus and S. tigrinus based on concatenated sequences (1892 bp) of the mitochondrial DNA COI and COII regions. Numbers at nodes indicate bootstrap support. (B) Maximum likelihood tree with the molecular clock assumption. Numerals indicate the estimated divergence times.
Biogeographic pattern of S.tigrinus
The phylogenetic tree of Fig. 6 shows that genetic differentiation of node 3 (lineage D + E and lineage F) and node 4 (lineage D and E) were estimated to have occurred during the glacial periods. In the glacial period during the Pleistocene, a number of populations became extinct, even as a number of haplotypes had vanished from populations due to decreasing of population size because of the colder climate on the Sea of Japan side than on the Pacific Ocean side. It may be inferred from the fact that the volume of inflow of the Tsushima Current known as a warm current was small during the glacial periods due to the low sea level (Kitamura and Kimoto, 2004).
However, a small number of populations with a few haplotypes in each lineage could survive in restricted areas such as refugia and represent the ancestral populations of lineages D, E and F, which had low genetic divergence (Table 4); in particular, lineages D and F are thought to have quite low population numbers. In the interglacial period, ancestral populations in lineages D and F began to rapidly expand population size and distribution range to the current size and distribution. This hypothesis was supported by a star-like haplotype networks and the significant negative values showed Tajima's D and Fu's Fs (Fig. 4 and Table 5). Therefore, the bottleneck effect seems to have an impact on lineages D and F. Moreover, lineage F has almost the same topology as the haplotype network for lineage D, being star-like with two or more step haplotypes that are more distant from the central haplotype (Hap_T31) than those of lineage D (Fig. 4), and fewer haplotypes were shared between populations than for lineage D that the range expansion of lineage F started earlier than for lineage D because of the warm climate and oceanic current (Tsushima current) in western Japan compared with cooler climate in northern Japan.
On the other hand, within lineage E, haplotypes were not shared between populations with one exception—Hap_ T12 occurred in the Ishikawa and Fukui populations (Fig. 4). Also, no major haplotypes occurred in this lineage. The haplotype network showed a few haplo-groups connected to each other by many missing haplotypes. Though these missing haplotypes might be not detected, these findings suggest the presence of stably maintained populations in each locality over a long time, few extinction events of populations, and rare events of gene flow between populations, particularly, between the Kyoto and Ishikawa + Fukui populations.
Comparison of genetic structures of S. griseus and S. tigrinus
The two closely related species, S. griseus and S. tigrinus, are similar to each other not only in morphology and habitat preference but also in the genetic characteristics in two ways: the division of three major lineages corresponding lineages A, B and C in S. griseus and lineages D, E and F in S. tigrinus (Figs. 3 and 4); and the largest proportion of the total genetic variance is among lineages, not among populations and not within populations (Table 7). Moreover, the two species also would have been simultaneously influenced in geographical history and the repeated climate changes with repeated glacial and interglacial periods during the Pleistocene. Nevertheless, those glacial and interglacial periods had different effects on the two species.
These two species are quite different from each other in their genetic structure S. Scepticus griseus showed relatively “high” genetic diversity (Tables 1 and 5) caused by little extinction of populations and repeated cycles of isolation of populations and gene flow between populations, resulting in a complicated genetic structure in this species. On the other hand, S. tigrinus showed relatively “low” genetic diversity (Tables 1 and 5) due to a number of population extinctions and the subsequent bottleneck effects, resulting in a more simple genetic structure in this species than in S. griseus. The identification of these two species, S. griseus and S. tigrinus, with different genetic structures provides a unique scenario for studying the evolutionary history in each of these species.
This analysis showed two distinct lineages between S. griseus and S. tigrinus inhabiting Pacific Ocean and Sea of Japan sides, respectively, corroborating the findings of Kudo et al. (2012) and Niikura et al. (2015). Gene flow between weevils on the Pacific Ocean and Sea of Japan sides rarely occurs. Niikura et al. (2015) also showed that the genetic structure of Tylos granuliferus, a semiterrestrial coastal isopod, on the main islands of Japan was strongly influenced by oceanic currents. However, this study of S. griseus and S. tigrinus revealed that the influences of geohistorical events on their genetic structure played a larger role than that of oceanic currents, as was observed for several other coastal species (Itani, 2000; Batch et al., 2004). In terms of biogeography, it is interesting to understand the genetic structure of other coastal species and to determine whether the genetic structure was influenced by oceanic currents or not.
In the present study of these coastal weevils, large genetic differences were observed within one species, and we also understood that gene flow might be rare between populations. Recently, some artificial sandy beaches were constructed along the coast due to increasing coastal erosion (Uda, 1997). As such, man-made alterations to the coastal regions may cause genetic pollution. Thus, it is important to understand the genetic structures of coastal species from a viewpoint of conservation biology in order to avoid genetic pollution.
ACKNOWLEDGMENTS
We thank Prof. S. Okajima (Tokyo University of Agriculture), Dr. T. Mita (Kyushu University) and Dr. H. Kato (The University of Tokyo) for their help in many ways. We are indebted for the generous offers of materials to the following former TUA members: Mr. Y. Kudo, Ms. S. Sato, Mr. Y. Ohno, Mr. S. Kojima, Mr. K. Kasai and Mr. T. Ichikawa. This study was supported in part by Grants-in-Aid from the JSPS KAKENHI (16K07484 to NK and 15K06937 to HK).
AUTHOR CONTRIBUTIONS
YY, H K and NK conceived and designed the study. YY, HK and TI collected samples. YY wrote the manuscript and performed the experiments. YY and TI prepared the Figures. NK and WD helped genetic analyses and revised this manuscript.
