Population Genetic Structure of the Bumblebee, <italic>Bombus ignitus</italic> (Hymenoptera: Apidae), Based on Mitochondrial COI Gene and Nuclear Ribosomal ITS2 Sequences
Population Genetic Structure of the Bumblebee, Bombus ignitus (Hymenoptera: Apidae), Based on Mitochondrial COI Gene and Nuclear Ribosomal ITS2 Sequences
International Journal of Industrial Entomology. 2013. Sep, 27(1): 142-158
Copyright ©2013, Korean Society of Sericultural Science
  • Received : July 16, 2013
  • Accepted : September 09, 2013
  • Published : September 30, 2013
Export by style
Cited by
About the Authors
Hyung Keun Oh
College of Agriculture & Life Sciences, Chonnam National University, Gwangju 500-757, Korea
Hyung Joo Yoon
Department of Agricultural Biology, National Academy of Agricultural Science, Suwon 441-100, Korea
Joo Young Lee
College of Agriculture & Life Sciences, Chonnam National University, Gwangju 500-757, Korea
Jeong Sun Park
College of Agriculture & Life Sciences, Chonnam National University, Gwangju 500-757, Korea
Iksoo Kim
College of Agriculture & Life Sciences, Chonnam National University, Gwangju 500-757, Korea

The bumblebee, Bombus ignitus (Hymenoptera: Apidae), is a valuable natural resource that is widely utilized for greenhouse pollination in South Korea. Understanding the magnitude of genetic diversity and geographic relationships is of fundamental importance for long term preservation and utilization. As a first step, we sequenced a partial COI gene of mitochondrial DNA (mtDNA) corresponding to the “DNA barcode” region and the complete internal transcribed spacer 2 (ITS2) of nuclear ribosomal DNA from 88 individuals collected in nine South Korean localities. The complete ITS2 sequences were longest among known insects, ranging in size from 2,034 bp ~ 2,052 bp, harboring two duplicated 112-bp long repeats. The 658-bp long mtDNA sequences provided only six haplotypes with a maximum sequence divergence of 0.61% (4 bp), whereas the ITS sequences provided 84 sequence types with a maximum sequence divergence of 1.02% (21 sites). The combination of the current COI data with those of published data suggest that the B. ignitus in South Korea and China are genetically a large group, but those in Japan can be roughly separated into another group. Overall, a very high per generation migration ratio, a very low level of genetic fixation, and no discernable hierarchical population were found to exist among the South Korean populations of B. ignitus , which suggests panmixia. This finding is consistent with our understanding of the dispersal capability of the species.
The European bumblebee, Bombus terrestris , has been commercially introduced in several parts of the world including South Korea and neighboring countries (Iwasaki, 1995; Mitsuhata, 2000). Despite its importance for crop pollination, the side effects, such as competition among bumblebee species (Inoue et al. , 2008) and disturbance of pollination of other bumblebees adapted to a specific crop (Dohzono et al. , 2008) have also been reported. In South Korea, an interspecific hybridization, genetic contamination, and competitive copulation between B. terrestris and B. ignitus have also been reported under laboratory conditions (Yoon et al. , 2009). Thus, a substantial effort has been directed at establishing an artificial mass rearing system for B. ignitus in order to substitute the B. terrestris (Yoon et al. , 2002).
Currently, B. ignitus is under commercialized by several companies in South Korea. Once field-released, the bumblebee populations that have been adapted to local environment may possibly be replaced by artificial ones, which could reduce local genetic diversity and fitness (Kawecki and Ebert, 2004). Therefore, knowledge on genetic diversity and geographic relationships of B. ignitus is essential for long-term preservation and artificial selection. Nevertheless, studies on bumblebee in this regard have been extremely limited for South Korean populations (Lee et al. , 2006). Under such circumstances, Tokoro et al. (2010) recently reported that B. ignitus occurring in Japan is slightly, but obviously different from those of South Korea and China by sequencing 1,048 bp of the CO gene of mitochondrial DNA (mtDNA). However, this conclusion was made on the basis of a single population from each Korea and China, and, thus, further confirmation may be required using additional populations.
In this study, therefore, we collected a total of 88 individuals of B. ignitus from nine South Korean localities, and sequenced 658 bp of the mitochondrial COI gene and complete internal transcribed spacer 2 (ITS2) from nuclear ribosomal RNA. The sequence information was subjected to the analysis of population genetic structure with the inference of genetic diversity and genetic relationships for South Korean populations. Further, COI sequences were combined with those of Tokoro et al. (2010) to better understand the geographic variation of the Korean populations relative to Japanese populations. The COI region we employed corresponded to the “DNA Barcode” region (Hebert et al. , 2003), which has been used to provide insight into the patterning of withinspecific genomic diversity (Hajibabaei et al. , 2007). The nuclear ITS2, which is located between the 5.8S and 28S rRNA genes, has been broadly investigated in population genetic and contemporary gene flow studies (Mukabayire et al. , 1999;
Lager Image
Sampling location of Bombus ignitus in Korea. Generallocality names and number of sequence indicated are as follows: 1,Jeongseon, Baekjeon-ri, Gangwondo Province; 2, Youngheungdo,Incheon City; 3, Chuncheon City, Gangwondo Province; 4,Busan City; 5, Suwon City, Gyeonggido Province; 6, Namyangju,Gyeonggido Province; 7, Taebaek City, Gangwondo Province;8, Jungseon, Gangwondo Province; and 9. Muju, JeollabukdoProvince.
Marcilla et al. , 2001).
Materials and Methods
- Insects
A total of 88 B. ignitus was collected from nine localities around the beautiful corydalis Colidalis speciosa (Fumariaceae) and the cherry blossom Prunus yedoensis (Rosaceae) in South Korea. Although most sampling was performed for only one day, those from Jeongseon (locality 1) were mingled with two days samples, resulting in the highest sampling size ( Table 1 ). Regardless of similar sampling effort, more individuals were collected successfully in some localities when compared with others. Sampling locality, number of individuals, date of collection, and GenBank accession numbers for individual COI and ITS2 sequences are provided in Table 1 and the locality map is shown in Fig. 1 .
Inventory of samples and summary of sequencing resultsWithin-parentheses indicate number of individuals collected at each locality and collection date.
Lager Image
Inventory of samples and summary of sequencing results Within-parentheses indicate number of individuals collected at each locality and collection date.
- DNA extraction, primer, PCR, and sequencing
Total DNA was extracted with a Wizard Genomic DNA Purification Kit, in accordance with the manufacturer's instructions (Promega, USA). A 658-bp region of the mitochondrial COI gene corresponding to the “DNA Barcode” region (Herbert et al. , 2003) was amplified using a pair of primer sets reported by Kim et al. (2009). For the amplification of ITS2, the primers NG02955 and AB052895 located on the 5.8S and 28S rDNAs, respectively, were successfully used (Ji et al. , 2003). In order to complete the entire ITS2 the internal primer BITSF, 5’-CGTAGTGTTGTCCTCGTGACCGA-3’, was designed from the conserved region of the B. ignitus ITS2 sequences.
For amplification of the COI gene, PCR was conducted under the following conditions: an initial denaturation step at 94℃ for 7 min, a 35-cycle amplification (94℃ for 1 min, 54~56℃ for 1 min, and 72℃ for 1 min), and the final extension step for 7 min at 72℃. For amplification of ITS2, PCR was conducted under the following conditions: an initial denaturation step at 94℃ for 7 min, a 35-cycle amplification (94℃ for 40 sec, 55~64℃ for 20 sec, and 72℃ for 20 sec), and the final extension step for 2 min at 72℃. Although COI amplicons were directly sequenced, those of ITS2 were cloned into a pGEM-T Easy vector (Promega, USA). XL1-Blue competent cells (Stratagene, USA) were transformed with the ligated DNA, and the resultant plasmid DNA from one clone per individual was isolated with a Plasmid Miniprep Kit (Dyne Bio Inc., Korea). Sequencing was performed using the ABI PRISM ® BigDye ® Terminator v3.1 Cycle Sequencing Kit with an ABI 3100 Genetic Analyzer (PE Applied Biosystems, USA).
- Sequence analysis
COI gene and ITS2 sequenc es were delimitated and aligned using MAFFT ver. 6 (Katoh et al. , 2002). The boundary decision to remove 5.8S at the 3′ end and 28S at the 5′ end of the ITS2 sequence was made using the Hidden Markov Modelbased ITS2 annotation software (Keller et al. ,, 2009). When the homologous sequences from two individuals differed by one or more nucleotide bases (for both COI and ITS2) or one insertion/deletion (indel) position (for ITS2), the sequences were considered different haplotypes (for COI) or sequence types (for ITS2). Haplotype or sequence type designations were applied to new sequences as they were discovered (i.e., BARBI01, BARBI02, and BARBI03 for COI, and ITSBI01, ITSBI02, ITSBI03 and so on for ITS2).
For the alignment of the indel-containing ITS2 region sequences, GBlocks version 0.91b software (Castresana, 2000) was used to select conserved regions for the subsequent phylogenetic and population level analyses. Resultantly, a total of 84 ITS2 sequence types with 2,034~2,052 positions, including gaps from 84 ITS2 sequence types with ~2,068 positions were obtained, conserving 96% of the original sequences.
- Phylogenetic analysis
To determine the relationships among the COI haplotypes, among ITS2 sequence types, and to detect any describable groups, phylogenetic analysis was conducted via the maximum-parsimony (MP) method (Fitch, 1971) using PAUP* (Phylogenetic Analysis Using Parsimony and Other Method*) ver. 4.0b10 software (Swofford, 2002) and Bayesian Inference (BI) method using MrBayes ver. 3.1 (Huelsenbeck and Ronquist, 2001), respectively. The analysis for the MP method was conducted by heuristic search. The reliability of the trees was assessed by 1,000 iterations of bootstrapping (Felsenstein, 1985).
For the BI analysis, the substitution model was selected by comparing Akaike information criterion scores (Akaike, 1974) using Modeltest ver. 3.7 (Posada and Crandall, 1998). The best-fit model selected was GTR (Lanave et al. , 1984) + I + G for the COI gene and K80 + G (Kimura, 1980) for ITS2. The confidence values of the BI tree were presented as the Bayesian posterior probabilities in percentages (BPP). The homologous regions of the within-generic species B. ardens were adapted, respectively, from Kim et al. (2009) for the COI gene and Oh et al. (2009) for ITS2 as an outgroup to root the trees.
- Network construction
Haplotype (COI) or sequence type (ITS2) relationships were determined using the median-joining algorithm (Bandelt et al. , 1999) incorporated in SplitsTree ver. 4.11.3 (Huson and Bryant, 2006). This method adds to the network median vectors (consensus sequences) by starting with the minimum spanning trees combined within a single network. Such vectors can be interpreted as possibly extinct unsampled sequences or extinct ancestral sequences (Bandelt et al. , 1999).
- Population genetic estimates
Haplotype and nucleotide diversity, both of which are reflective of genetic diversity within each locality, were determined using Arlequin ver. 3.5 (Excoffier and Lischer, 2010), whereas maximum sequence divergence within each locality was estimated by extracting within-locality estimates of unrooted pairwise distances from PAUP* ver. 4.0b (Swofford, 2002). The population pairwise genetic distance ( F ST ) and a permutation test for its significance (1,000 bootstraps) were evaluated using Arlequin ver. 3.5 (Excoffier and Lischer, 2010). Pairwise Nm (the product of the effective population size, Ne , and migration rate, m ) values were employed to estimate the pairwise F ST based on the equilibrium relationship: F ST = 1 / (2 Nm + 1) for the COI gene and F ST = 1 / (4 Nm + 1) for ITS2. For these estimates only populations with more than two individuals, possessing more than two haplotypes for the COI gene or more than two sequence types for ITS2 were subjected to analyses. Genetic relationships among populations and sets of populations were assessed by the Holsinger and Mason-Gamer (H-MG) method (1996). A detailed rationale of this method is described in the original study of Holsinger and Mason-Gamer (1996).
- Sequence analysis
A total of six COI haplotypes (BARBI01 ~ BARBI06) and 84 ITS2 sequence types (ITSBI01 ~ ITSBI84) were obtained by sequencing 88 individuals of B. ignitus ( Table 1 ). Although the COI gene contained no indels, presenting all identical 658 bp, ITS2 was variable in length, ranging from 2,034 to 2,052 bp. The length of B. ignitus ITS2 was the longest among known insects as far as available data are considered ( Table 2 ). Such a long ITS2 sequence in B. ignitus derives, in part, from the presence of a total of 112-bp long, two repeat units, repeated at the beginning region within ~420 bp of B. ignitus ITS2 ( Fig. 2 ). An uncorrected pairwise comparison between pairs of haplotypes demonstrated that the sequence divergence ranged from 0.15 to 0.61% (1 ~ 4 bp) among the six haplotypes ( Table 3 ). From the ITS2 sequences a total of 84 sequence types, a sequence divergence ranging from 0.04% to 1.02% was obtained (one ~21 positions; data not shown). In contrast to the higher sequence type diversity of ITS2, the maximum sequence divergence among sequence types was only slightly higher compared to the COI gene (1.02% vs. 0.61%).
Among the six haplotypes, four were found in a single locality as a single individual (BIBAR03, BIBAR04, BIBAR05, and BIBAR06), but the haplotype BIBAR02 was found in four localities (localities 1, 2, 7, and 8) in a total of five individuals and BIBAR01 was found in seven localities in a total of 79 individuals, accounting for 89.8% of the samples ( Table 1 ). Thus, the distribution of B. ignitus haplotypes can be summarized as a restricted local distribution in most haplotypes, with a wide distribution only in a limited number of haplotypes. On the other hand, no ITS2 sequence type was found in more than one
Summery of ITS2 size and G+C contents in other insects
Lager Image
Summery of ITS2 size and G+C contents in other insects
Lager Image
Beginning sequences of Bumbus ignitus ITS2 rDNA, which contained two repeat units.
Pairwise comparisons among six haplotypes obtained from mitochondrial COI gene sequence ofBombus ignitusNumbers above the diagonal are mean distance values; numbers below the diagonal are absolute distance values.
Lager Image
Pairwise comparisons among six haplotypes obtained from mitochondrial COI gene sequence of Bombus ignitus Numbers above the diagonal are mean distance values; numbers below the diagonal are absolute distance values.
Lager Image
Phylogenetic trees obtained from six mitochondrial COIhaplotypes of the Bombus ignitus. (A) Tree acquired via the MPmethod incorporated in the PAUP ver. 4.0b10 software (Swofford,2002). The numbers on the branches represent bootstrap values of1,000 replications. (B) Tree obtained via the Bayesian Inferencephylogram. Numbers at each node specify the BPP. Withinparenthesesindicate the locality number, from which the particularsequence type was obtained. Bombus ardens was used as anoutgroup in order to root tree. The scale bar indicates the number ofsubstitutions per site.
locality, except for the sequence type ITSBI01, which was found in five localities as a single individual at each locality (1, 2, 5, 7, and 8) ( Table 1 ).
- Phylogenetic relationships
Phylogenetic analysis was conducted to determine the relationships among the six haplotypes of B. ignitus obtained in this study ( Fig. 3 ). All haplotypes were weakly associated or unresolved, owing principally to low genetic divergence among them. One exception was the clustering between BIBAR01 and BIBAR04 in the BI analysis with slightly higher node support (BPP = 0.82), but this group was very weakly supported in the MP analysis as 55%.
Previously, Tokoro et al. (2010) sequenced 1,048-bp of the COI gene from B. ignitus collected from a Korean and Chinese locality and several Japanese localities. They obtained a total of 15 haplotypes, which consisted of six obtained from both Korea and China and nine from Japan and found that Japanese haplotypes were different from both South Korean and Chinese haplotypes using the MP and Neighbor-Joining method. In order to prepare a total dataset for phylogenetic analysis, the current six haplotypes obtained in this study were combined and aligned with those of Tokoro et al. (2010), resulting in 13 new haplotypes with 463 bp in length. These included five from this study (BIBAR01, BIBAR02 = BIBAR03, BIBAR04, BIBAR05, and BIBAR06) and eight from Tokoro et al. (2010) (F1 = F2 = F3, F4, F5 = F6, J1 = J2 = J3 = J5, J4, J6, J7 = J8, J9). Phylogenetic analysis by MP and BI methods demonstrated that the haplotypes found in South Korea and China formed one relatively inclusive group (0.92 by BI and 52% by MP), whereas those found in Japan did not form a monophyletic group (data not shown). This
Within-locality diversity estimates ofBombus ignitusobtain from COI and ITS2 sequencea) Sample sizeb) Number of haplotypesc) Haplotype diversityd) Number of polymorphic sitese) Maximum sequence divergencef) Mean number of pairwise differencesg) Nucleotide diversity
Lager Image
Within-locality diversity estimates of Bombus ignitus obtain from COI and ITS2 sequence a) Sample size b) Number of haplotypes c) Haplotype diversity d) Number of polymorphic sites e) Maximum sequence divergence f) Mean number of pairwise differences g) Nucleotide diversity
Lager Image
Median-Joining networks indicating the relationship amongCOI haplotypes obtained in this study and those of Tokoro et al.(2010). Haplotype in group A were found in Korean and Chineselocalities, whereas those in group B were found solely in Japan.Identical haplotypes generated from the selection of overlapping463 nucleotides between this study and that of Tokoro et al. (2010)are indicated next to the haplotype names within parentheses.Locality numbers and names are presented under the haplotypetype names within parentheses. The branch lengths representthe amount of character-state changes occurring on that branch.Circular dots represent haplotypes found in this study, whereasrectangles indicate the hypothetical ones that were not found inthis study.
was likely due to the short sequence length (463 bp) employed in this study.
In order to further understand the relationships among B. ignitus haplotypes, the COI-based median-joining network was prepared using the haplotypes obtained in this study and those of Tokoro et al. (2010). Among the five haplotypes found in this study (BIBAR01, BIBAR02 = BIBAR03, BIBAR04, BIBAR05, and BIBAR06) the haplotypes BIBAR01 and BIBAR04 were somewhat distantly related to others, but these were grouped together with other haplotypes found in the South Korea and China with some complexity, due to the presence of more than one most parsimonious connection among haplotypes, forming one distinct group (termed group A) ( Fig. 4 ). On the other hand, the haplotypes found solely in Japan (J1 = J2 = J3 = J5, J4, J6, J7 = J8, J9) formed a somewhat distant independent group (termed group B) that was derived from a single founder, J1 ( Fig. 4 ).
Phylogenetic analyses among the 88 ITS2 rDNA sequence types of B. ignitus have shown largely unresolved multibranches, but a few sequence types were grouped together, although the nodal support for the grouping was relatively low ( Fig. 5 ). MP-based analysis has shown that the majority of the sequence types were weakly associated or unresolved, but two groups were exceptional, with bootstrap values of either 85% or 99% ( Fig. 5 A). Similarly, BI-based analysis also resulted in largely unresolved or weakly associated groups ( Fig. 5 B). In an attempt to further understand the relationships among B. ignitus sequence types, an ITS2-based median-joining network was prepared ( Fig. 6 ). In contrast to the MP- and BI-based
Lager Image
Phylogenetic trees obtained from 84ITS2 sequence types of the Bombus ignitus. (A) Tree acquired via the MP method incorporatedin the PAUP ver. 4.0b10 software (Swofford, 2002). The numbers on the branches represent bootstrap values of 1,000 replications. (B) Treeobtained via the Bayesian Inference phylogram. Numbers at each node specify the BPP. Within-parentheses indicate the locality number,from which the particular sequence type was obtained. Bombus ardens was used as an outgroup in order to root tree. The scale bar indicatesthe number of substitutions per site.
phylogenies, the network resulted in a larger number of inclusive groups, providing several star phylogenies, wherein sequence types were composed of centered ones and their derived ones. In any star phylogeny, the derived sequence types were composed of several sequence types that were originated from a diverse locality. Thus, the clustering pattern of the sequence types derived from these localities may reflect non-locality-based clustering, indicating gene flow.
- Population genetic estimates
Within-locality diversity was estimated in terms of haplotype
Lager Image
Median-Joining networks indicating the relationship among 84 ITS2 sequence types. The branch lengths represent the amount ofcharacter-state changes occurring on that branch. Circular dots represent sequence types found in this study, whereas rectangles indicate thehypothetical ones that were not found in this study.
diversity ( H ), maximum sequence divergence (MSD), mean number of pairwise differences (MPD), and nucleotide diversity (π) for both COI and ITS2 sequences from the localities with a sample size and haplotype number of ≥ two. Consequently, the three and five localities for the COI and ITS2 sequences, respectively, were subjected to analysis ( Table 4 ). In the case of the COI gene, all three localities displayed substantially low estimates of π, ranging from 0.000767 (locality 1; Jeongseon) to 0.000217 (locality 8; Muju). Considering the maximum difference between the highest and lowest estimates only accounted for about 3.5-fold, no single particular locality could be attributed to a center for genetic diversity. Similarly, the ITS2 sequences also provided only 1.75-fold of the maximum difference between the highest and lowest estimates of π, ranging from 0.004615 (locality 8; Muju) to 0.002645 (locality 2; Youngheungdo). Thus, the ITS2 data also indicated that no single particular locality could be attributed to a center for genetic diversity either. On the other hand, the values of π in the ITS2 were about five to 20-fold higher than those of the COI gene sequence.
Both COI and ITS2-based genetic distance ( F ST ) revealed
Fixation indices (FST) and migration rate (Nm) between pairs of populations ofBombus ignitusobtain from mitochondrial COI gene*p < 0.05.inf, infinite.
Lager Image
Fixation indices (FST) and migration rate (Nm) between pairs of populations of Bombus ignitus obtain from mitochondrial COI gene * p < 0.05. inf, infinite.
no statistically significant F ST in any case (Online Resources 9 and 10). Consistent with the F ST data, the estimate of pergeneration migration rates ( Nm ) showed an overall high gene flow between pairs of populations both in the COI and ITS2 sequences ( Tables 5 and 6 ). To better understand the nature of the genetic relationships between the B. ignitus populations, the hierarchical relationships among localities were analyzed
Fixation indices (FST) and migration rate (Nm) between pairs of populations ofBombus ignitusobtain from ITS2 sequence*p < 0.05.inf, infinite.
Lager Image
Fixation indices (FST) and migration rate (Nm) between pairs of populations of Bombus ignitus obtain from ITS2 sequence * p < 0.05. inf, infinite.
Lager Image
Hierarchical relationships among localities analyzed usingthe Holsinger and Mason-Gamer method (1996). The dendrogramobtained from (A) COI and (B) ITS2. The value at each node isthe distance (D) between its two daughter nodes and the p valueis the significance of differentiation (based on 10,000 randomresamplings).
( Fig. 7 ). The COI-based analysis demonstrated that there was a close relationship between Jeongseon (locality 1) and Youngheungdo (locality 2), although this relationship was not statistically significant ( p > 0.05). In the case of the ITS2-based analysis, Busan and Muju (localities 4 and 8) formed into one group, and Jeongseon and Ulleungdo (localities 1 and 9) formed into another group, and the Busan + Muju group clustered together with the Youngheungdo (locality 2), forming (Youngheungdo + (Busan + Muju)). Nevertheless, separation into these groups was not statistically supported ( p > 0.05). Thus, hierarchical analyses also indicated nongeographic distance-based population clustering.
- Major characteristic of each molecule
Other available bee species present in South Korea have a COI sequence divergence of 0.3% for the within-generic bumblebee B. ardens (Kim et al. , 2009) and 0.76% for the mason bee Osmia cornifrons (Kim et al. , 2008), indicating that the sequence divergence was < 1%. Thus, the maximum sequence divergence of the COI gene found as 0.61% in B. ignitus was moderate when compared with other bee species occurring in South Korea.
To the best of our knowledge, the B. ignitus ITS2 sequence is the longest ever reported in insects, which ranged between 2,034 ~ 2,052 bp. Previously, two within-familial species, Melipona , were reported to harbor a ITS2 with lengths of 1,728 bp and 1,789 bp, respectively (De La Rúa et al. , 2007) and the length of the within-generic species, B. ardens ITS2 sequence was reported to be 1,971 ~ 1,984 bp (Oh et al. , 2009), but our B. ignitus was at least 50 bp longer than those of the B. ardens ITS2 sequences. Except for the Apidae family, the length of the ITS2 in other hymenopteran and other ordinal groups were far less than 1,000 bp in most cases ( Table 2 ). Similarly to B. ignitus , the B. ardens ITS2 sequences was 114-bp in length and contained four repeat units, repeated twice at the beginning region within ~430 bp of ITS2, but no such repeat sequence was detected in the ITS2 regions of any other insect, including two within-familial Melipona species (De La Rúa et al. , 2007). Thus, this repeat sequence of the Bombus ITS2 seems to be a molecular synapomorphic characteristic for the genus, although more Bombus species should be analyzed.
The sequence type diversity of ITS2 was also noteworthy in that nearly all individuals accounted for a different sequence type ( Table 1 ). One of the possible reasons for such high diversity in sequence type is a higher evolutionary ratio of ITS2 since this ITS2 may have been relatively free of structural and functional constrains (Tang et al. , 1996).
- Relationships between South Korean and Japanese B. ignitus
Previously, Shao et al. (2004) sequenced 450 bp of mitochondrial cytochrome b gene from 16 individuals collected from three localities in China, South Korea, and Japan, and no sequence variation was observed among them. However, further variable microsatellite analysis from 124 individuals showed statistically significant genetic differentiation of the Chinese and South Korean populations from the Japanese populations, which was possibly due to the geographical isolation of Japan (Shao et al. , 2004). Tokoro et al. (2010) also found that the genetics of Japanese populations were differentiated from South Korean and Chinese populations when a 1,048 bp of the COI gene was sequenced. Our median-joining network using the data from Tokoro et al. (2010) and the current study also suggested that there was a genetic differentiation between mainland populations and Japan ( Fig. 7 ). Collectively, this data indicate that the B. ignitus in South Korea and China are genetically a large group, but those in Japan is roughly separable into another group. This result reinforces the finding of Tokoro et al. (2010), even though limited sequence information was utilized. It is likely that the geographic connection between the Korean peninsula and mainland China facilitated gene glow, but the geographically isolated Japanese populations may have had little or no chance for gene flow during continental separation at about 110,000 ~ 280,000 ( Kawamura, 2006 ).
- Genetic relationships among South Korean populations
Our F ST and Nm analyses on the basis of both COI and ITS2 indicate that B. ignitus populations occurring in South Korea were not genetically differentiated from each other ( Tables 5 and 6 ). Furthermore, the hierarchical relationship data based on both the COI and ITS2 indicate that all B. ignitus populations in South Korea were genetically interrelated ( Fig. 7 ). Consistent with the population genetic perspective, the phylogenetic analyses ( Figs. 3 and 5 ) and the median-joining networks ( Figs. 4 and 6 ) among South Korean populations also indicated the importance of gene flow in the B. ignitus .
Although no precise report on the dispersal distance of B. ignitus is available it has been known that the within-generic species, B. terrestris , B. lucorum , and B. lapidaries are capable of dispersal for several thousand kilometers (Mikkola, 1978). Furthermore, it has been theorized that a small number of migrants may be sufficient to prevent drift for genetically homogeneous populations (Hartl and Clark, 1989). Considering these findings, the B. ignitus populations in South Korea seem to be a single panmictic population, regardless of geographic distance.
Summarized, we collected a total of 88 individuals of B. ignitus from nine South Korean localities and sequenced 658 bp of the mitochondrial COI gene and 2,034 ~ 2,052-bp long complete ITS2 rDNA. The sequence analysis of B. ignitus COI provided a moderate to low magnitude of sequence divergence, whereas that of ITS2 indicated higher variability, presenting a total of 84 sequence types. In particular, the lengths of ITS2 were the longest among known insects. Considering the F ST and Nm estimates, along with phylogenetic analyses, it appears that the geographic populations of B. ignitus have a high ratio of gene flow, which indicates the potential/actual dispersal ability of the species without long-term zoogeographic barriers. This was likely due to high dispersal ability.
This work was supported by a grant (Code 20070401-034-004) from the BioGreen 21 Program, Rural DevelopmentAdministration, Republic of Korea.
Akaike H (1974) A new look at the statistical model identification. IEEE Trans Autom Contr 19 716 - 723    DOI : 10.1109/TAC.1974.1100705
Alam MT , Bora H , Das MK , Sharma YD (2008) The type andmysorensis forms of the Anopheles stephensi (Diptera: Culicidae) inIndia exhibit identical ribosomal DNA ITS2 and domain-3 sequences. Parasitol Res 103 75 - 80    DOI : 10.1007/s00436-008-0930-7
Bandelt HJ , Forster P , Röhl A (1999) Median-joining networks forinferring intraspecific phylogenies. Mol Biol Evol 6 37 - 48    DOI : 10.1093/oxfordjournals.molbev.a026036
Beebe NW , Ellis JT , Cooper RD , Saul A (1999) DNA sequence analysisof the ribosomal DNA ITS2 region for the Anopheles punctulatusgroup of mosquitoes. Insect Mol Biol 8 381 - 390    DOI : 10.1046/j.1365-2583.1999.83127.x
Castresana J (2000) Selection of conserved blocks from multiple alignments for their use in phylogenetic tool. Curr Opin Genet Dev 8 669 - 674
De La Rúa P , De May-Itzá JW , Serrano J , Quezada-Euán JJG (2007) Sequence and RFLP analysis of the ITS2 ribosomal DNA in twoNeotropical social bees, Melipona beecheii and Melipona yucatanica(Apidae, Meliponini). Insectes Sociaux 54 418 - 423    DOI : 10.1007/s00040-007-0962-5
Dohzono I , Kunitake YK , Yokoyama J , Goka K (2008) Alienbumblebee affects native plant reproduction through interactions withnative bumblebee. Ecology 89 3082 - 3092    DOI : 10.1890/07-1491.1
Erasmus JC , Van Noort S , Jousselin E , Greeff JM (2007) Molecularphylogeny of fig wasp pollinators (Agaonidae, Hymenoptera) ofFicus section Galoglychia. Zool Scripta 36 61 - 78    DOI : 10.1111/j.1463-6409.2007.00259.x
Excoffier L , Lischer HEL (2010) Arlequin suite ver 3.5: a new series ofprograms to perform population genetics analyses under Linux andWindows. Mol Ecol Resour 10 564 - 567    DOI : 10.1111/j.1755-0998.2010.02847.x
Felsenstein J (1985) Confidence limits on phylogenics: an approachusing the bootstrap. Evolution 29 783 - 791    DOI : 10.2307/2408678
Fitch WM (1971) Toward defining the course of evolution: minimalchange for a specific tree topology. Syst Zool 20 406 - 416    DOI : 10.2307/2412116
Gómez-Zurita J , Juan C , Petitpierre E (2000) Sequence, secondarystructure and phylogenetic analyses of the ribosomal internaltranscribed spacer 2 (ITS2) in the Timarcha leaf beetles (Coleoptera:Chrysomelidae). Insect Mol Biol 9 591 - 604    DOI : 10.1046/j.1365-2583.2000.00223.x
Hackett BJ , Gimnig J , Guelbeogo W , Constantini C , Koekemoer LL (2000) Ribosomal DNA internal transcribed spacer (ITS2) sequencesdifferentiate Anopheles funestus and Anopheles and A. rivulorum ,and uncover a cryptic taxon. Insect Mol Biol 9 369 - 374    DOI : 10.1046/j.1365-2583.2000.00198.x
Hajibabaei M , Singer GA , Hebert PD , Hickey DA (2007) DNAbarcoding: how it complements taxonomy, molecular phylogenetics and population genetics. Trends Genet 23 167 - 172    DOI : 10.1016/j.tig.2007.02.001
Hartl DL , Clark AG (1989) Principles of population genetics. Sinauer Sunderland, Mass
Hebert PDN , Cywinska A , Ball SL , DeWaard JR (2003) Biologicalidentifications through DNA barcodes. Proc R Soc Lond B Biol Sci 7 313 - 321    DOI : 10.1098/rspb.2002.2218
Holsinger KE , Mason-Gamer RJ (1996) Hierarchical analysis ofnucleotide diversity in geographically structured populations. Genetics 142 629 - 639
Huelsenbeck JP , Ronquist F (2001) MrBayes: Bayesian inference ofphylogeny. Bioinformatics 17 754 - 755    DOI : 10.1093/bioinformatics/17.8.754
Hung YT , Chen CA , Wu WJ , Shih CJ (2004) Phylogenetic utility ofthe ribosomal internal transcribed spacer 2 in Strumigenys spp.(Hymenoptera: Formicidae). Mol Phylogenet Evol 32 407 - 415    DOI : 10.1016/j.ympev.2004.03.010
Huson DH , Bryant D (2006) Application of phylogenetic networks inevolutionary studies. Mol Biol Evol 23 254 - 267    DOI : 10.1093/molbev/msj030
Inoue MN , Yokoyama J , Washitani I (2008) Displacement of Japanesenative bumblebee by the recently introduced Bombus terrestris (L.)(Hymenoptera, Apidae). J Insect Conserv 12 135 - 146    DOI : 10.1007/s10841-007-9071-z
Iwasaki M (1995) Introduction of commercial bumblebees into Japan. Honeybee Sci 16 17 - 21
Ji Y , Zhang D , He L (2003) Evolutionary conservation and versatility ofa new set of primers for amplifying the ribosomal internal transcribedspacer regions in insects and other invertebrates. Mol Ecol Notes 3 581 - 585    DOI : 10.1046/j.1471-8286.2003.00519.x
Katoh K , Misawa K , Kuma K , Miyata T (2002) MAFFT: a novelmethod for rapid multiple sequence alignment based on fast Fouriertransform. Nucleic Acids Res 30 3059 - 3066    DOI : 10.1093/nar/gkf436
Kawamura Y (2006) Recent progress in paleontological studies on theQuaternary mammals of Japan. Mammal Sci 47 107 - 114
Kawecki TJ , Ebert P (2004) Conceptual issues in local adaptation. Ecol Lett 7 1225 - 1241    DOI : 10.1111/j.1461-0248.2004.00684.x
Keller A , Schleicher T , Schultz J , Muller T , Dandekar T , Wolf M (2009) 5.8S-28S rRNA interaction and HMM-based ITS2 annotation. Gene 430 50 - 57    DOI : 10.1016/j.gene.2008.10.012
Kim HY , Lee KY , Lee SB , Kim SR , Hong MY , Kim DY , Kim I (2008) Mitochondrial DNA sequence variation of the Mason bee, Osmiacornifrons (Hymenoptera: Apidae). Int J Indust Entomol 16 75 - 86
Kim MJ , Yoon HJ , Im HH , Jeong HU , Kim MI , Kim SR , Kim I (2009) Mitochondrial DNA sequence variation of the bumblebee, Bumbus ardens (Hymenoptera: Apidae). J Asia-Pacific Entomol 12 133 - 139    DOI : 10.1016/j.aspen.2009.02.003
Kimura M (1980) A sample method for estimating evolutionary rateof base substitution through comparative studies of nucleotide sequences. J Mol Evol 16 111 - 120    DOI : 10.1007/BF01731581
Lanave C , Preparata G , Saccone C , Serio G (1984) A new method forcalculating evolutionary substitution rates. J Mol Evol 20 86 - 93    DOI : 10.1007/BF02101990
Lee ML , Yoon HJ , Jin BR (2006) Morphological and DNA variabilitiesof Bombus ignitus Smith (Apodae : Hymenoptera). Kor J Apicul 21 139 - 144    DOI : 10.3348/kjr.2006.7.2.139
Marcilla A , Bargues MD , Ramsey JM , Magallon-Gastelum E , Salazar-Schettino PM , Abad-Franch F , Dujardin JP , Schofield CJ , Mas-Coma S (2001) The ITS-2 of the nuclear rDNA as a molecular marker forpopulations, species and phylogenetic relationships in Triatominae(Hemiptera: Reduviidae), vectors of Chagas disease. Mol Phylogenet Evol 18 136 - 142    DOI : 10.1006/mpev.2000.0864
Marinucci M , Romi R , Mancini P , Di Luca M , Severini C (1999) Phylogenetic relationships of seven palearctic members of themaculipennis complex inferred from ITS2 sequence analysis. Insect Mol Biol 8 469 - 480    DOI : 10.1046/j.1365-2583.1999.00140.x
Mikkola K (1978) Spring migrations of wasps and bumble bees onsouthern coast of Finland (Hymenoptera Vespidae and Apidae). Ann Entomol Fennici 44 10 - 26
Mitsuhata M (2000) Pollination of crops with bumblebee colonies inJapan. Honeybee Sci 21 17 - 25
Muccio T , Marinucci M , Frusteri L , Maroli M , Pesson B , Gramiccia M (2000) Phylogenetic analysis of Phlebotomus species belongingto the subgenus Larrousius (Diptera: Psychodidae) by ITS2 rDNAsequences. Insect Biochem Mol Biol 30 387 - 393    DOI : 10.1016/S0965-1748(00)00012-6
Mukabayire O , Boccolini D , Lochouarn L , Fontenille D , Besansky NJ (1999) Mitochondrial and ribosomal internal transcribed spacer (ITS2)diversity of the African malaria vector Anopheles funestus. Mol Ecol 8 289 - 297    DOI : 10.1046/j.1365-294X.1999.00567.x
Oh HK , Yoon HJ , Kim MJ , Jeong HU , Kim SR , Hwang JS , Bae CH , Kim I (2009) ITS2 ribosomal DNA sequence variation ofthe bumblebee, Bombus ardens (Hymenoptera: Apidae). Genes Genomics 31 293 - 303    DOI : 10.1007/BF03191202
Posada D , Crandall KA (1998) Modeltest: testing the model of DNAsubstitution. Bioinformatics 14 817 - 818    DOI : 10.1093/bioinformatics/14.9.817
Samara R , Monje JC , Reineke A , Zebit CPW (2008) Genetic divergenceof Trichogramma aurosum Sugonjaev and Sorokina (Hymenoptera:Trichogrammatidae) individuals based on ITS2 and AFLP analysis. J Appl Entomol 132 230 - 238    DOI : 10.1111/j.1439-0418.2007.01245.x
Shao ZY , Mao HX , Fu WJ , Ono M , Wang DS , Bonizzoni M , Zhang YP (2004) Genetic structure of Asian population of Bombus ignitus. J Hered 95 46 - 52    DOI : 10.1093/jhered/esh008
Stouthamer R , Gai Y , Koopmanschap AB , Platner GR , Pinto JD (2000) ITS-2 sequences do not differ for the closely related speciesTrichogramma minutum and T. platneri. Entomol Exp Appl 95 105 - 111    DOI : 10.1046/j.1570-7458.2000.00647.x
Swofford DL (2002) PAUP*. Phylogenetic Analysis Using Parsimony(*and other method) software ver. 4.10 Sinauer Associates Sunderland, Massachusetts
Tang JM , Toe L , Black C , Unnasch TR (1996) Intraspecic heterogenietyof rDNA internal transcribed spacer in the Simu-lium damnosum(Diptera: Simulidae) complex. Mol Biol Evol 13 244 - 252    DOI : 10.1093/oxfordjournals.molbev.a025561
Thanwisai A , Kuvangkadilok C , Baimai V (2006) Molecular phylogenyof black flies (Diptera: Simuliidae) from Thailand, using ITS2 rDNA. Genetica 128 177 - 204    DOI : 10.1007/s10709-005-5702-z
Tokoro S , Yoneda M , Kunitake YK , Goka K (2010) Geographicvariation in mitochondrial DNA of Bombus ignitus (Hymenoptera,Apidae). Appl Entomol Zool 45 77 - 87    DOI : 10.1303/aez.2010.77
Vobis M , Haese JD , Mehlhorn H , Mencke N , Blagburn BL , Bond R , Denholm I , Dryden MW , Payne P , Rust MK , Schroeder I , Vaughn MB , Bledsoe D (2004) Molecular phylogeny of isolates ofCtenocephalides felis and related species based on analysis of ITS1,ITS2 and mitochondrial 16S rDNA sequences and random bindingprimers. Parasitol Res 94 219 - 226    DOI : 10.1007/s00436-004-1201-x
Wagener B , Reineke A , Loehr B , Zebitz CPW (2006) Phylogeneticstudy of Diadegma species (Hymenoptera: Ichneumonidae) inferredfrom analysis of mitochondrial and nuclear DNA sequences. Biol Control 37 131 - 140    DOI : 10.1016/j.biocontrol.2006.01.004
Weekers PHH , De Jonckheere JF , Dumont HJ (2001) Phylogenetic relationships inferred from ribosomal ITS sequences and biogeographic patterns in representatives of the genus Calopteryx (Insecta: Odonata) of the West Mediterranean and adjacent West European zone. Mol Phylogenet Evol 20 89 - 99    DOI : 10.1006/mpev.2001.0947
Yoon HJ , Kim SY , Lee KY , Lee SB , Park IG , Kim I (2009) Interspecific hybridization of the bumblebees Bombus ignitus and B. terrestris. Int J Indust Entomol 18 41 - 48
Yoon HJ , Kim SE , Kim YS (2002) Temperature and humidity favorable for colony development of the indoor-reared bumblebee, Bombus ignitus. Appl Entomol 37 419 - 423    DOI : 10.1303/aez.2002.419