INTRODUCTION
Parasitic lice (Insecta: Psocodea: ‘Phthiraptera’) are obligate parasites of mammals and birds and spend their entire lifecycle on the hosts. Because of this permanent association with their hosts, lice have been used as a model system for studies of host-parasite co-divergence (Page, Reference Page2002). Several studies comparing louse and host phylogenies have provided evidence for co-divergence between lice and their hosts (e.g., Hafner et al. Reference Hafner, Sudman, Villablanca, Spradling, Demastes and Nadler1994; Page et al. Reference Page, Lee, Becher, Griffiths and Clayton1998; Johnson et al. Reference Johnson, Williams, Drown, Adams and Clayton2002b; Light and Hafner, Reference Light and Hafner2008; Smith et al. Reference Smith, Light and Durden2008). In contrast, several cases in which substantial incongruence between host and parasite trees exists have been identified (Johnson et al. Reference Johnson, Adams and Clayton2002a, Reference Johnson, Williams, Drown, Adams and Claytonb; Weckstein, Reference Weckstein2004; Banks et al. Reference Banks, Palma and Paterson2006). All of these studies have focused on species-level phylogenies, and there are only a few studies focusing on co-divergence within species, among populations (Whiteman et al. Reference Whiteman, Kimball and Parker2007; Štefka et al. Reference Štefka, Hoeck, Keller and Smith2011). Thus, it is still generally unclear whether co-divergence at the species level is driven by population level co-divergence that precedes speciation. Population genetic and phylogeographic studies will be an important bridge between microevolution and co-phylogenetic patterns (Johnson et al. Reference Johnson, Adams, Page and Clayton2003a).
To understand the causes of congruence and incongruence between host and parasite phylogenies, careful consideration of host and parasite behaviours, effects of other members of the parasite community (e.g. phoresis by hippoboscid flies), and geographical information are needed (Johnson et al. Reference Johnson, Adams and Clayton2002a, Reference Johnson, Williams, Drown, Adams and Claytonb). As long as these factors are accounted for, parasite phylogeny could also provide insights regarding the host's phylogeny (Hopkins, Reference Hopkins1942; Reed et al. Reference Reed, Smith, Rogers, Hammond and Clayton2004; Johnson et al. Reference Johnson, Kennedy and McCracken2006; Whiteman et al. Reference Whiteman, Kimball and Parker2007). In particular, the molecular substitution rate for lice is generally much faster than that of their hosts (Hafner et al. Reference Hafner, Sudman, Villablanca, Spradling, Demastes and Nadler1994; Huelsenbeck et al. Reference Huelsenbeck, Rannala and Yang1997; Page, Reference Page1996; Page et al. Reference Page, Lee, Becher, Griffiths and Clayton1998; Johnson et al. Reference Johnson, Cruickshank, Adams, Smith, Page and Clayton2003b), and substitutional differences may accumulate between diverging louse populations at shallow time-scales even though their hosts have not yet accumulated any mutational differences. Therefore, population genetic and phylogeographic studies of parasitic lice also have the potential to uncover hidden population structure in their hosts.
In addition to studies of co-divergence, lice have also been used as an important system in studying morphological evolution of co-evolving systems. The evolution of morphological traits can be affected by many factors, such as phylogenetic constraints, gene flow, genetic drift, sexual selection, and, most importantly, selection pressure from the environment. Generally, it is very difficult to separate these factors clearly. Parasitic lice spend their entire life cycle on the body of the host and thus the environmental context for their morphological evolution can be more easily identified than that for free-living organisms. For this reason, lice are good models for studying the impact of the host on morphological evolution of these parasites. Previous studies have shown interesting co-evolutionary patterns between lice and their hosts. For example, Reed et al. (Reference Reed, Hafner and Allen2000) clearly showed that the head groove width of gopher chewing lice is strongly correlated with the hair diameter of the host. Mammal lice hang onto the host hair using the head groove, so that a mismatch between groove size and hair diameter likely increases the risk to the louse of falling off or being groomed off the host. Johnson et al. (Reference Johnson, Bush and Clayton2005) also found that the body width of wing lice from pigeons and doves is strongly correlated with host body size. Host body size dictates the width of the feather interbarb space, in which wing lice insert themselves to escape host-preening defences. As in the case of the gopher lice, a mismatch between the body width of the louse and the feather interbarb space increases the risk of being preened off the host (Clayton et al. Reference Clayton, Bush, Goates and Johnson2003; Bush and Clayton, Reference Bush and Clayton2006).
Most previous studies of louse morphological evolution have focused only on functionally significant characters, such as louse body width vs feather interbarb space or louse head groove size vs hair diameter. However, the functional significance (if any) of many taxonomically and phylogenetically useful morphological characters in lice is not clear. In addition, recent progress on the higher level molecular phylogeny of parasitic lice has revealed that morphological convergence can be very frequent for characters used in louse taxonomy and phylogeny reconstruction (Johnson et al. Reference Johnson, Yoshizawa and Smith2004, Reference Johnson, Shreve and Smith2012; Smith et al. Reference Smith, Page and Johnson2004; Yoshizawa and Johnson, Reference Yoshizawa and Johnson2006, Reference Yoshizawa and Johnson2010). For lice, little is known about how these taxonomically and phylogenetically useful morphological characters vary across populations and how these characters evolve over micro-evolutionary time-scales. Such knowledge would be useful to understand the evolution of morphology over macro-evolutionary time. Quantitative comparisons of morphometric data and phylogeny are an effective way to study these characters (Klingenberg and Gidaszewski, Reference Klingenberg and Gidaszewski2010), but no such study has been done for lice to date.
Here, we undertake a co-phylogeographic study of sika deer and their parasitic lice. The sika deer (Cervus nippon) is distributed in East Asia and throughout the Japanese Archipelago. Japanese populations are subdivided into 6 subspecies and are highly variable in body size (from 30 to 120 kg in male body mass: Ohtaishi, Reference Ohtaishi1986). Such a large variation in body size within a species of mammal is exceptionally rare (Ohdachi et al. Reference Ohdachi, Ishibashi, Iwasa and Saito2009). However, morphology-based subspecific subdivisions within sika deer are not supported by mitochondrial phylogeny, and 2 distinct mitochondrial lineages have been identified for Japanese sika deer. Sika deer host a parthenogenetic chewing louse species, Damalinia sika, making this a potentially ideal system to evaluate patterns of co-phylogeography and morphological evolution. As mentioned above, morphological evolution can be affected by many factors, including gene flow via interbreeding. Damalinia sika is an obligately parthenogenetic louse, and thus gene flow via interbreeding among populations can be excluded as a factor affecting morphological evolution, making this species useful for distinguishing other factors impacting morphological evolution of this species.
In this study, we estimated the mitochondrial population level phylogenies of sika deer and their lice, based on sequences of the COI gene. Co-phylogeographic patterns were investigated by comparing tree topologies and genetic divergence between hosts and parasites. We also discuss novel insights regarding deer population structure based on phylogeographic patterns observed in their lice. To examine patterns of louse morphological evolution, head shape was analysed using morphometric analyses. We compared the louse population tree from COI sequences with the results of cluster analysis of morphometric data.
MATERIALS AND METHODS
Samples
Deer and deer louse samples were collected throughout the Japanese archipelago (Fig. 1). Most louse samples were collected by the authors and fixed in 99% ethanol immediately. Hairs of deer were also sampled with lice and were stored in 99% ethanol. Some deer samples were collected by local hunters, and frozen skins were sent to the authors. Lice were collected from the skins and fixed with 99% ethanol. Numbers of deer and lice samples examined are listed in Table 1. In most cases, samples of lice came from the individual deer included in this study. However, the Taiwanese population of louse was represented from a host individual that was unavailable for sequencing. Sequences from the host in this case (C. n. taiounus) was represented by sequences obtained from GenBank. In addition, sequences from the Far East Russian population of sika deer were obtained from GenBank and included in the phylogeographic analysis of deer, even though louse samples were not available for analysis. This was done to stabilize the tree search by breaking the long branch leading to Taiwanese sika deer (Table 1). Cervus elaphus and Damalinia forficula were used as outgroups for the host and parasite trees, respectively. These species are not associated, but are simply used for rooting the ingroup topologies because they are among the closest relatives of C. nippon and D. sika.
Fig. 1. Map of collecting sites for sika deer and lice. Circles indicate the collection sites, with different colours showing different louse clades. For numbers of deer and lice samples from each locality, see Table 1.
Table 1. List of the samples used in this study
Sequences of a portion of the mitochondrial COI gene were obtained for both deer (838 bp) and lice (776 bp). Though the lengths of these sequences were slightly different, the sequence for the deer completely overlapped that for the lice, making these regions homologous. Therefore, the sequence data can be utilized for the comparison of host-parasite substitution rates (see below). The mitochondrial control region is known to evolve rapidly and thus is commonly used for phylogeography of mammals (Nagata et al. Reference Nagata, Masuda, Tamate, Hamasaki, Ochiai, Asada, Tatsuzawa, Suda, Tado and Yoshida1999). However, because the mitochondrial genome of deer lice is composed of minicircles (Cameron et al. Reference Cameron, Yoshizawa, Mizukoshi, Whiting and Johnson2011), homology of the control regions between deer and deer lice cannot be assumed. No nuclear markers for phylogeographic studies are known for deer (Barbosa and Carranza, Reference Barbosa and Carranza2010). Therefore, to compare homologous gene regions between hosts and parasites, COI is the best option to date. In addition, because sika deer lice are parthenogenetic, mitochondrial DNA is one of the best markers for tracking the history of these parthenogenetic lines. Procedures for total DNA extraction, PCR, and sequencing followed those described by Yoshizawa and Johnson (Reference Yoshizawa and Johnson2010), but primers were newly designed here (see Supplementary data, online version only). Alignment was done by eye because of the lack of indels in this protein-coding gene. The aligned matrix in nexus format and csv file for statistical analyses (see below) are available as Supplementary data, (online version only). Sequence data also are deposited at GenBank (Table 1). Uncorrected and likelihood-based genetic distances were calculated using PAUP* 4b10 (Swofford, Reference Swofford2002).
Phylogeographic analyses
Phylogeographic analyses were performed with the parsimony (PAUP*), likelihood (PAUP* and PhyML 3.0: Guindon et al. Reference Guindon, Dufayard, Lefort, Anisimova, Hordijk and Gascuel2010), and Bayesian (MrBayes 3.12: Ronquist and Huelsenbeck, Reference Ronquist and Huelsenbeck2003) criteria. Best-fit models for the ML and Bayesian analyses were calculated by hLRT using Modeltest (Posada and Crandall, Reference Posada and Crandall1998) and MrModeltest (Nylander, Reference Nylander2004), respectively (Supplementary data, online version only). For ML analyses, using PAUP* a heuristic search with TBR branch swapping was implemented for each data set, using a NJ starting tree. For MP analysis, 100 random addition replicates using TBR heuristic branch swapping were performed. PAUP* was also used for MP bootstrapping, with 1000 replicates. PhyML was used to conduct likelihood-based 1000 bootstrap replicates, using NNI branch swapping, each with Bio-NJ as a starting tree. For Bayesian analyses, we performed 2 runs, each with 4 chains for 10 000 000 generations, and trees were sampled every 1000 generations. For each run, the first 5000 trees were excluded as a burnin and the 50% majority consensus tree of the remaining trees was calculated to estimate posterior probabilities. We also estimated ParsimonySplits networks as implemented in SplitsTree4 (Huson and Bryant, Reference Huson and Bryant2006) using default settings. The outgroups were excluded from the splits analyses.
Co-divergence analyses
Co-divergence analysis was performed using the Jungle algorithm (Charleston, Reference Charleston1998) in TreeMap 2.02b (Page and Charleston, Reference Page and Charleston1998). Samples with identical haplotypes were trimmed to a single terminal taxon. All non-co-divergence events were given equal costs of 1 and co-divergence was given a cost of 0. As discussed below, a number of non-co-divergence events were identified from the Jungle analysis. To test whether there was significant evidence for a lack of strict co-divergence, some alternative genealogical topologies were tested using the approximately unbiased test (AU test: Shimodaira, Reference Shimodaira2002) using CONSEL (Shimodaira and Hasegaea, Reference Shimodaira and Hasegawa2001). The partition homogeneity test (Farris et al. Reference Farris, Kallersjo, Kluge and Bult1994, Reference Farris, Kallesjo, Kluge and Bult1995) was also performed using PAUP* to evaluate whether host and parasite data sets represented significantly different topologies (Johnson et al. Reference Johnson, Drown and Clayton2001). To perform the test, a data matrix in which hosts and parasites were matched in 1:1 associations (as described by Johnson et al. Reference Johnson, Drown and Clayton2001) was prepared by excluding all the samples lacking a host-parasite association (i.e., outgroups and sika deer from Ussuri) and also eliminating the samples for which both the hosts and parasites from an identical locality show exactly identical haplotypes. Using this 1:1 data set, 1000 replicates of the partition homogeneity tests were performed each with TBR branch swapping. This 1:1 association data set was also used to compare the genetic distance between deer and lice. PAUP* was used to estimate the genetic distances, and uncorrected pairwise distances of total data, uncorrected pairwise distances of 4-fold degenerate sites, and likelihood-based distances using the models estimated above were calculated. MEGA5 (Tamura et al. Reference Tamura, Peterson, Peterson, Stecher, Nei and Kumar2011) was used to identify 4-fold degenerate sites.
Analyses of louse morphology
For parasites, host-imposed selection may be more important than the parasite population or phylogenetic history for shaping the patterns of morphological evolution. To evaluate the correspondence between louse morphology and population structure of both deer and lice, the shape of the louse head was analysed. Elliptic Fourier analysis (Kuhl and Giardiana, Reference Kuhl and Giardina1982) as implemented in the software package SHAPE (Iwata and Ukai, Reference Iwata and Ukai2002) was used. Morphological terminology follows Lyal (Reference Lyal1985).
For morphological analysis, 138 louse specimens from 15 localities were examined, and each locality was treated as a terminal. The morphological samples did not include those used for DNA extraction during which Proteinase K caused breakage of the head capsule. Specimens that were used only for morphological analyses were cleared using 10% KOH solution for 1 night at room temperature. All specimens were rinsed with 80% ethanol (3 times), dehydrated with 99·5% ethanol (2 times), and slide mounted with Euparol.
The head capsule was photographed through the microscope using a digital camera. Using the digital image, the head outline was traced and scanned. Because the head shape of D. sika is close to a circle, the initial alignment of head shape had to be re-orientated. To do this, the right margin of the outline was oriented on the right margin of ventral carina. Following the method of Smith (Reference Smith2000), 20 harmonics were calculated from the digitized outline, while excluding effects of size. The resulting 80 coefficients were subjected to principal component analysis using SHAPE. The 13 principal component scores (PC 1–13) that represented over 95% of the variation were used for canonical discriminant analysis. We then calculated Mahalanobis distances between each group of specimens from the same locality. Then, using the Mahalanobis distances among locality groups, we performed cluster analysis using the group average method. Discriminant and cluster analyses were performed using JMP 9.0 (SAS, 2011).
RESULTS
Phylogeography of deer and deer lice
ML analysis of each data set recovered a single most likely tree (Figs 2 and 3). MP analyses recovered 18 (for deer) and 11 (for louse) trees of equal length, and one of these was the same as the tree estimated from the ML analyses. ParsimonySplits networks are provided as Supplementary Figs 1 and 2, (online version only).
Fig. 2. Phylogenetic tree of sika deer individuals estimated using maximum likelihood. Branch lengths are proportional to substitutions per site. Numbers associated with each branch indicate MP bootstrap/ ML bootstrap/ Bayesian posterior probabilities. Shades indicate haplotype group names used for Figs 4 and 5. Numbers in circles indicate co-divergence points corresponding to those in Fig. 3. See also Figs 4–6.
Fig. 3. Phylogenetic tree of deer lice estimated using maximum likelihood. Branch lengths are proportional to substitutions per site. Numbers associated with each branch indicate MP bootstrap/ ML bootstrap/ Bayesian posterior probabilities. Shades indicate haplotype group names used for Figs 4 and 5. Numbers in circles indicate co-divergence points corresponding to those in Fig. 2. Broken-line circles indicate samples from Kyushu and adjunct islands. See also Figs 4–6.
Three major lineages (North, South, and Continental + Taiwanese Clades) were recovered for sika deer (Fig. 2, Suppl. Fig. 1, online version only). Support values for all clades were generally strong (>77% BS for all clades and 98% PP for South Clade), but only low Bayesian support values were obtained for North (98–100% BS vs 82% PP) and Continental + Taiwanese Clades (77–83 BS vs 81% PP). Two Japanese lineages (Northern and Southern lineages) composed a clade, and the Continental + Taiwanese samples formed the other clade. However, resolution of the basal relationships among these lineages was weakly supported (59% BS and 56% PP). Low resolution of the basal divergences was also evident from the splits network (Suppl. Fig. 1, online version only). Analysis of the smaller data set (1:1 data, see Materials and Methods) placed Taiwanese sika deer as the sister of the Southern clade, making Japanese sika deer paraphyletic (Supplementary data, online version only). The Southern Clade included all samples from Kyushu and its adjacent islands. The Northern Clade consisted of the samples from Hokkaido and Honshu.
The analyses of the louse data set recovered 2 distinct lineages (Clades I and II: Fig. 3 and Suppl. Fig. 2, online version only). Support values for both clades were very strong. Clade I was subdivided into 2 clades (Ia and Ib). Clade Ia was well supported and included samples from Hokkaido, northern Honshu and Fukuoka. The basal branch of Clade Ib was extremely short and weakly supported. This clade was composed of samples from Nikko, southern Honshu, Tsushima and Kagoshima. In the split network, samples of Tsushima were placed closer to Clade Ia than the other samples of Clade Ib, although the network connecting them was extremely short (Suppl. Fig. 2, online version only). Clade II was composed of the samples from southern Kyushu and adjacent islands, and the sample from Taiwan.
Co-divergence analysis
Initial Jungle analysis using TreeMap 2 recovered 4 equally optimal results with 4 co-divergence events and 13 non-co-divergence events (Figs 4 and 5). A standard randomization test showed that the number of co-divergence events was significantly higher than that expected by chance in 1000 randomizations of the parasite tree (P = 0·005). However, most co-phylogenetic analyses do not accommodate parasites that occur on multiple host terminals, such as a single species of parasites on multiple host species or parasite populations found on multiple host populations (Johnson et al. Reference Johnson, Adams and Clayton2002a). Analytically these are often post hoc considered to be recent host switching events. Because of these difficulties in dealing with widespread parasites, the results from the TreeMap analysis could not account for 2 shallow host switching events of lice from deer in Fukuoka to North1 and to Toya populations (asterisks in Fig. 5). There still was significant congruence between host and parasite trees (P = 0·034) even if these two additional non-co-divergence events were allowed.
Fig. 4. Tanglegram of deer and louse trees with lines connecting host-parasite associations.
Fig. 5. Optimal Jungle reconstruction using TreeMap2. Note that two host switch events indicated by asterisks were added by hand, because of the inability of Jungles analysis to deal with parasite lineages occurring on multiple host lineages. Labels for louse are distinguished from deer by a “p” heading. See text for detailed explanation.
Because louse phylogeography did not closely correspond to host deer divergence patterns, we also examined these topological incongruences statistically. When the distinction between the South and North populations in the deer was forced as a constraint on the louse tree (excluding the Taiwanese deer louse, which was kept within the Southern Clade, together with Kyushu deer louse, as in the original louse phylogeny), the resulting tree was significantly worse than the unconstrained tree using the AU test (P < 0·001). In contrast, when only the basal divergence between the Taiwanese and Japanese populations as evidenced in the deer phylogeny was constrained in the louse tree, this constraint tree was not significantly worse than the original using the AU test (P = 0·217). Similarly for the deer, when the two major clades identified by the phylogenetic analyses of the louse was constrained in the tree for sika deer (excluding the Taiwanese sika deer, which was kept separate from the Japanese population as in the original deer phylogeny), the constrained tree was significantly worse using the AU test (P = 0·027). Likewise, the partition homogeneity test recovered significant heterogeneity between the host-parasite COI sequences (P = 0·001).
Comparisons of genetic distances
Fig. 6 shows plots of the uncorrected p-distances between pairs of lice and those from their deer hosts. These plots contain roughly 5 clusters of points (Groups A–E). The divergence events of deer and lice corresponding to the plot groups are indicated by shade in the figure, and the plots corresponding to the co-divergence events are indicated by circled numbers (Figs 2, 3, 5 and 6). The uncorrected p-distances of 4-fold degenerate sites and likelihood-based distances also showed identical trends (not shown). When points in groups B and C were excluded (because they clearly represent some incongruence between louse and deer population histories), and correlation between deer and louse genetic distances were calculated, the substitution rate of the COI gene for lice was estimated to be about 4·8 times faster than that of deer (Fig. 6). This estimate also accords with the fact that the highest divergences among lice are around 12%, while the highest among deer are only 2·5%. Although plot group E corresponds to the co-divergence event between Kumamoto and Yakushima + Kuchinoerabu populations (Fig. 5), these points are also outliers of the estimated slope (Fig. 6).
Fig. 6. Plots of uncorrected pairwise distances estimated from deer and louse COI sequences. Shadings indicate deer and lice divergence events as noted on left and lower margins. Dashed outlines and associated letters indicate plot groups as discussed in the text. Numbers in circles indicate co-divergence points corresponding to those in Figs 2, 3 and 5.
Morphological variation of lice
Except for the populations from Kinkazan and Nara, all clusters identified by morphological analysis corresponded to the lineages recovered from the louse molecular phylogenetic analysis (Figs 7 and 8). Note that molecular phylogenetic analysis divided the Nikko population into 2 distinct lineages (Fig. 3), but the majority of these specimens were placed in the Clade Ia, which corresponded to the cluster identified by morphological analysis. In contrast, except for 2 populations collected from the same subspecies of deer (i.e., C. n. yakushimae: note that lice from these two populations are also genetically closely related: Fig. 3), the clusters identified from louse head morphology did not show any grouping corresponding to deer subspecies or geography.
Fig. 7. Illustration of shape variation in louse head morphology, with variation contained within each of the first four principal components indicated.
Fig. 8. Dendrogram resulting from cluster analysis of the louse head shape principal components. Terminal labels indicate louse collection site (left), their host subspecies (represented by first three letters of subspecific epithet: middle), and the mitochondrial lineages for each louse population (right).
PC1 explained variation in the shape of conus, depth of the antennal groove, head width, and width of the pulvinus (Fig. 7: proportion = 32·6%). PC2 also explained variation in the protrusion of the antennal groove margin, protrusion of the temple, shape of the conus, and width of the pulvinus (proportion = 26·9%). PC3 explained variation in the antennal groove (proportion = 10·1%). PC4 explained clear variation in the depth of the pulvinus, and explained some variation in the depth of the antennal groove (proportion = 9·1%). The proportion of morphological variation explained by PC5-13 was less than 5% (not shown). However, PC5 (proportion = 4·7%) still contained detectable variation in the depth of the antennal groove, and PC6 (proportion = 2·7%) explained clear variation in the length of the conus. PC7-13 explained very little of the overall variation.
DISCUSSION
Co-divergence between sika deer and lice
Comparison of deer and louse trees using Jungle analyses recovered significant patterns of topological congruence. However, these analyses only account for the branching pattern and do not take into account the relative timing of branching events. Using information from the genetic divergence plots, 3 of 4 reconstructed co-divergent events are also concordant with the relative timing of branching events (1, 3 and 4 in Fig. 6), providing reciprocal support for the results of both types of analyses. However, plot Group E, which would correspond to a co-divergence between Kumamoto and Yakushima + Kuchinoerabu population (2 in Fig. 6), is apparently an outlier for genetic distance as compared to the other co-divergence events. Therefore, even though the branching patterns of these populations might indicate co-divergence, based on genetic distances, this is more likely a recent host switch. When one of the co-divergence events was excluded (i.e., number of co-divergence events was set to 3 and one more host switching event was allowed), the tree randomization test showed that the number of co-divergence events was not significantly higher than that expected by chance (P = 0·19). Significant mismatches between the deer and louse trees are also evident from the plots of the genetic distances, the AU test, and the partition homogeneity test.
In contrast, the AU test did not reject the possibility of basal diversification of the Taiwanese deer louse from the Japanese population. Alternatively, when the smaller data set was analysed, the Taiwanese sika deer was placed to the sister of the Southern Clade, agreeing with the louse tree as far as the placement of Taiwanese population is concerned. However, there is a significant discordance between pairwise genetic distance comparisons for these deer and louse populations (plot group C in Fig. 6). Therefore, by considering both topology and genetic distances, the close relationship between louse populations of Taiwan and southern Kyushu is unlikely to be a result of their co-divergence and, as estimated from Jungle analyses, this should also be interpreted as a result of a recent host switch.
The louse populations in Kyushu are highly structured genetically, even within Clade I, whereas those of Hokkaido and Honshu show much less genetic variability. Jungle analysis suggests that the louse populations in Hokkaido and Honshu originated via a recent host switch from deer in Kyushu. The deer populations in Honshu and Hokkaido show some genetic structure, but the genetic variation of lice in Honshu and Hokkaido is apparently less structured and less variable, even though the substitution rates of lice are much faster than those for deer. The lack of genetic variability and structure in these louse populations also supports the hypothesis, from Jungles analysis, that they originated from a recent host switch, which may have been associated with a genetic bottleneck.
The evidence for multiple host switching events suggests that deer lice migrate among deer populations quite readily and independent of deer population structure. Mismatches between host and louse phylogenies at species level are quite frequent for bird lice (Johnson et al. Reference Johnson, Adams and Clayton2002a, Reference Johnson, Williams, Drown, Adams and Claytonb; Weckstein, Reference Weckstein2004; Banks et al. Reference Banks, Palma and Paterson2006). In contrast, all previous studies of mammal lice have identified highly significant co-divergence patterns between mammal lice and their hosts at the species (Hafner et al. Reference Hafner, Sudman, Villablanca, Spradling, Demastes and Nadler1994; Light and Hafner, Reference Light and Hafner2007; Light et al. Reference Light and Hafner2008; Smith et al. Reference Smith, Light and Durden2008) or population levels (Demastes et al. Reference Demastes, Spradling, Hafner, Spies, Hafner and Light2012). Although trends similar to other mammal lice were expected in these deer lice, it appears that different dynamics may operate for the deer-lice system, and these are explored below.
Phylogeography of deer lice
The most basal divergence within deer lice was estimated to correspond to the divergence between continental and Japanese deer lineages. Clade II appears to have its ancestor in Taiwanese populations, and Japanese populations seem to have originated by a recent host switch from the Taiwanese population to the South Clade. The common ancestor of Clade I was estimated to parasitize the South deer clade, and the louse populations on the North deer clade seem to have originated by a very recent host switch from louse populations parasitizing the South deer clade.
This scenario is reasonable based on the additional evidence provided by genetic distances. Deer lice collected from Kyushu contain haplotypes from all major louse clades (Clade Ia, Ib and II) and, even within Clade I, they are quite variable genetically (shown by dotted circles: Fig. 3). In contrast, although lice from Hokkaido and Honshu are also divided into 2 subclades (Clade Ia and Ib), their haplotypes are either identical with (Clade Ia) or only slightly different from (Clade Ib) those of Kyushu. Therefore, recent host switches of deer lice from Kyushu to Hokkaido + Honshu populations are highly plausible.
This scenario implies 2 possibilities. The first possibility is that the common ancestor of the Northern deer clade might not have hosted any louse. The Northern deer clade is thought to have colonized the Japanese Archipelago during the last ice age via Sakhalin (see below). If deer lice are not tolerant of cold temperatures, this scenario might be plausible. Although no quantitative data are available, we only collected 2 lice from 2 deer at Nishiokoppe, the northern most locality of the present sampling area, and the fieldwork was conducted in mid-winter. Based on other collections we have made, this was an exceptionally small number of lice infecting these deer. A second possibility is that the louse populations originally parasitizing deer in Hokkaido and Honshu were completely replaced by those from Kyushu following the recent host switch events. Recently, Damalinia sika was introduced to the USA and has since become widespread and is replacing native lice very rapidly in some native deer populations (Mertins et al. Reference Mertins, Mortenson, Bernatowicz and Briggs Hall2011, unpublished data). This may indicate that certain lineages of deer lice, particularly parthenogenetic species such as Damalinia sika, may have higher population growth rates, giving them the potential to replace other populations.
Implications of louse population structure for deer phylogeography
Previous studies from the mitochondrial D-loop and CytB sequences show that Japanese sika deer possess 2 genetically distinct lineages (Tamate et al. Reference Tamate, Tatsuzawa, Suda, Izawa, Doi, Sunagawa, Miyahara and Tado1998; Nagata et al. Reference Nagata, Masuda, Tamate, Hamasaki, Ochiai, Asada, Tatsuzawa, Suda, Tado and Yoshida1999). Using sequences of the COI gene, the present analysis also supported this previous result. However, as discussed above, louse phylogeography is highly incongruent with deer population structure, even though one would expect that dispersal of these permanent parasites is highly tied to that of their hosts.
One possible explanation for the topological incongruence may be phoresis of these lice on hippoboscid flies as has been convincingly demonstrated for pigeon lice (Harbison et al. Reference Harbison, Jacobsen and Clayton2009; Harbison and Clayton, Reference Harbison and Clayton2011). Unlike lice, hippoboscid flies have wings and can disperse between host individuals relatively freely. Sika deer host 3 hippobiscid flies, Lipoptena fortisetosa, L. sikae and L. japonica (Mogi, Reference Mogi1975; Mogi et al. Reference Mogi, Mano and Sawada2002). The final instar larvae of these species leave the deer for pupation and, after emergence, louse-free adult flies fly to new hosts. Soon after the arrival to new host deer, the flies cut off their wings and thus are not likely to disperse to a new host (Yamauchi and Nakayama, Reference Yamauchi and Nakayama2006; Yamauchi, personal communication). In addition, there are no records of Damalinia sika attached to hippoboscid flies. Therefore, although the possibility cannot be excluded completely, phoresis does not seem a likely explanation for the topological incongruence between deer and louse trees.
The more likely explanation for incongruence between deer and louse trees is that lice migrate via male deer. The 2 geographically distinct lineages of Japanese sika deer have only been recovered using maternally inherited mitochondrial markers, which only represents the genetic structure of female deer (Nagata et al. Reference Nagata, Masuda, Tamate, Hamasaki, Ochiai, Asada, Tatsuzawa, Suda, Tado and Yoshida1999). In contrast, 2 such geographically distinct groups are not recovered in analyses of microsatellite markers (Goodman et al. Reference Goodman, Tamate, Wilson, Nagata, Tatsuzawa, Swanson, Pemberton and McCullough2001). In fact, the genetic structure of Japanese sika deer recovered using microsatellite markers is more congruent with the louse phylogeny, because the deer populations from Kyushu are divided into 2 clusters based on microsatellites. Furthermore, microsatellite data clusters a deer population from northern Kyushu together with populations from northern Honshu and Hokkaido (Goodman et al. Reference Goodman, Tamate, Wilson, Nagata, Tatsuzawa, Swanson, Pemberton and McCullough2001), which is also concordant with the louse phylogeny. Sika deer are polygamous, and females compose maternal groups that do not migrate and tend to stay within a male's territory. In contrast, young male deer (3–7 years old) migrate widely before obtaining a harem (Miura, Reference Miura1984; Takatsuki, Reference Takatsuki2006). Therefore, deer lice could be dispersed widely by male deer, and the dispersal of males is likely to be the main source for the incongruence between deer and lice mitochondrial population histories. Observations from deer nuclear markers suggest that the louse phylogeny does at least partially reflect the deer's population structure, with respect to the fact that male gene flow is quite frequent between 2 mitochondrial lineages of sika deer. The louse population structure also provides further support for the idea that subspecific division in sika deer does not reliably reflect population subdivision.
The 2 Japanese mitochondrial lineages of sika deer have been estimated to have already diverged on the continent before their invasion of the Japanese Archipelago (around 0·3–0·5 Mya), based on mitochondrial genetic distances and geographical history of the region (Nagata et al. Reference Nagata, Masuda, Tamate, Hamasaki, Ochiai, Asada, Tatsuzawa, Suda, Tado and Yoshida1999; Okada, Reference Okada, Takatsuki and Yamagiwa2008). The Southern Clade is hypothesized to have invaded via the Korean Peninsular in the mid-Pleistocene (−0·17 Mya), and the Northern Clade is estimated to have invaded via Sakhalin in late Pleistocene (60 000–10 000 ya). However, as mentioned above, these estimates are based only on the female line genetic structure. In the present analysis, the Taiwanese deer louse is placed very close to the deer louse from southern Kyushu, which implies that southern Japanese populations of sika deer may have had contact with continental deer more recently than estimated from the mitochondrial phylogeny of the deer alone.
As pointed out in previous studies (Johnson et al. Reference Johnson, Adams and Clayton2002a; Clayton and Johnson, Reference Clayton and Johnson2003; Weckstein, Reference Weckstein2004; Poulin et al. Reference Poulin, Krasnov, Mouillot and Thieltges2011), detailed information about louse and host behaviours is important for understanding the causes of similarities and discordance between the host and louse phylogenies. The behaviour of the sika deer is well studied (reviewed by Takatsuki, Reference Takatsuki2006), which greatly aided in interpreting the discordance between the deer and louse phylogenies. Alternatively, studies of louse phylogeography could also prove important for understanding the host's population structure. The substitution rate of louse mitochondrial genes is much faster than that of their vertebrate hosts (see below), and these insects can be transmitted both paternally as well as maternally. In the case of sika deer, good nuclear markers for population genetic studies, including those on Y-chromosome, have not been discovered to date (Barbosa and Carranza, Reference Barbosa and Carranza2010), and long distance dispersal of deer lice apart from their hosts is quite unlikely. In such cases, louse population genetic trees can be good indicators of hidden population structure in their hosts (Reed et al. Reference Reed, Smith, Rogers, Hammond and Clayton2004; Whiteman et al. Reference Whiteman, Kimball and Parker2007; Zohdy et al. Reference Zohdy, Kemp, Durden, Wright and Jernvall2012). This notion seems to be supported in our study by the fact that the louse tree is more similar to the deer population structure estimated from nuclear microsatellite markers (Goodman et al. Reference Goodman, Tamate, Wilson, Nagata, Tatsuzawa, Swanson, Pemberton and McCullough2001) than to that estimated from mitochondrial genes (Nagata et al. Reference Nagata, Masuda, Tamate, Hamasaki, Ochiai, Asada, Tatsuzawa, Suda, Tado and Yoshida1999).
Comparisons of substitution rates
By plotting pairwise uncorrected genetic distance comparisons for lice against those for their deer, roughly 5 groups of points are evident. Group A represents the shallowest divergences within both deer and lice. The phylogenetic trees do not show significant evidence for unevenness of substitution rates among populations for either deer or lice. Therefore, one of the slopes connecting A–B, A–E–C or A–D likely estimates the relative rate of substitution for lice vs deer of the COI gene. As discussed above, plot groups B and C clearly represent some incongruence between louse and deer population history. Given these considerations, the slope of the line connecting groups of comparisons A and D seems to be the best estimate for the relative rate of substitution between lice and deer (Fig. 6). Most of the comparisons between lineages involved in co-divergence (see below and Results) are also in Groups A and D.
This slope estimates the relative substitution rate to be around 4·8 times faster in lice than in deer when uncorrected pairwise distances are used. Comparisons of 4-fold degenerate sites (i.e. neutral sites: 137 bp for deer and 131 bp for lice) also provided a comparable value (4·9 times). This relative rate is greater than those estimated previously (2–4 times faster) for lice and their hosts (Hafner et al. Reference Hafner, Sudman, Villablanca, Spradling, Demastes and Nadler1994; Page, Reference Page1996; Huelsenbeck et al. Reference Huelsenbeck, Rannala and Yang1997; Page et al. Reference Page, Lee, Becher, Griffiths and Clayton1998; Reed et al. Reference Reed, Smith, Rogers, Hammond and Clayton2004; Light and Hafner, Reference Light and Hafner2007; Demastes et al. Reference Demastes, Spradling, Hafner, Spies, Hafner and Light2012). Phylogenetic studies show that mitochondrial substitution rate of deer lice is faster than that of other lice (Yoshizawa and Johnson, Reference Yoshizawa and Johnson2003, Reference Yoshizawa and Johnson2010: note that the species is mislabelled as Bovicola tibialis in Yoshizawa and Johnson, Reference Yoshizawa and Johnson2003: see Mertins et al. Reference Mertins, Mortenson, Bernatowicz and Briggs Hall2011 for further information). Therefore, the presently estimated relative rate (4·8 times) can be regarded as reasonable. Although a far greater relative rate (11 times) was estimated by comparing the branch lengths of likelihood trees that were congruent between doves and their wing lice (Johnson et al. Reference Johnson, Cruickshank, Adams, Smith, Page and Clayton2003b), a smaller estimate of relative rate of 1·6 times was obtained based on parsimony. Although not shown, comparisons of likelihood-based pairwise distances for deer and louse COI sequences resulted in a greater relative substitution rate (7·1 times). This value does not account for topology because of the many non-co-divergence events identified in the present analyses. Estimates based on branch lengths in a tree topology probably provide higher relative rates because more multiple substitutions can be accounted for using a tree-based approach, as was done for doves and their wing lice (Johnson et al. Reference Johnson, Cruickshank, Adams, Smith, Page and Clayton2003b). As mentioned above, the substitution rates of deer louse mitochondrial genes are faster than other lice sequenced to date (Yoshizawa and Johnson, Reference Yoshizawa and Johnson2003, Reference Yoshizawa and Johnson2010), possibly related to its parthenogenesis (Bengtsson, Reference Bengtsson2003; Henry et al. 2011; Shreve et al. Reference Shreve, Mockford and Johnson2011).
Micro-evolution of louse head morphology
The dendrogram obtained by cluster analysis of louse head morphology was highly concordant with the louse phylogenetic tree. The only obvious discordance was a difference in the placement of populations from Nara and Kinkazan, otherwise all the major clusters were concordant with the louse molecular tree. Thus, the shape of the louse head capsule appears to contain significant phylogenetic signal.
The morphology of the conus and neighbouring structures shows the greatest variation and these variables are the major factors in PCs 1–3 and major or minor factor in PCs 4–6. The function of the conus is poorly understood. However, this structure often shows sexual dimorphism and is frequently absent in male Trichodectidae (Lyal, Reference Lyal1985) so that the conus may not be a highly functional structure and thus selective pressure on this character may be relatively weak. The second structure showing large variation is the anterior margin of the pulvinus. This character is one of the factors in PCs 1, 2 and the major factor in PC 4. This structure is involved in grasping the host hair and thus is functionally important. Previous study has shown that this character is strongly correlated with host size and thus is likely to be under strong selection pressure (Reed et al. Reference Reed, Hafner and Allen2000).
The factors explaining the largest variation in louse head morphology appear to be both relatively neutral as well as characters expected to be under stronger selection. Variation in these characters clearly contains signal concordant with the louse phylogeographic history. In this case, the neutral characters appeared to comprise a larger fraction of the variation compared to those expected to be more functionally significant. This suggests that random drift of more neutral characters plays a major role in the micro-evolution of morphology for these lice. The evolution of functionally significant characters is slower compared to that of less functional characters, perhaps because of stabilizing selection. Interestingly, the shape of the conus is frequently used in taxonomic and phylogenetic studies, but the shape of the pulvinus is rarely used (e.g., Lyal, Reference Lyal1985; Smith, Reference Smith2001). A high level of variation and correlation with phylogeny of taxonomically important characters, such as dorsal head pattern and abdominal pigmentation, have also been identified in the louse genus Brueelia (Johnson et al. Reference Johnson, Adams and Clayton2002a), although the functional significance of these characters is less clear. Many important diagnostic characters used in louse taxonomy possibly originated by the accumulation of more neutral variation along independent phylogenetic lineages.
The present analysis shows that micro-evolutionary changes in deer louse morphology accumulate along phylogenetic lines. In contrast, Smith et al. (Reference Smith, Page and Johnson2004) and Johnson et al. (Reference Johnson, Shreve and Smith2012) showed that, across bird feather lice, convergence is frequent in functionally significant characters, such as body width (i.e., body vs wing lice) or the shape of the anterior head margin. Such convergence is evidence that natural selection produced the variation seen in louse morphology across bird lice. The morphology of deer lice also provides a potential example of convergence caused by similar selection pressures. Although most morphometric clusters are concordant with the louse phylogenetic tree, populations from Kinkazan and Nara were placed at discordant positions. Interestingly, sika deer of these populations have unusual histories. Kinkazan is a small island (9·6 km2), and the population size of sika deer on this island has been highly variable. Given this population fluctuation in their hosts, it is also likely that the lice would have experienced population bottlenecks. Dwarfishness of deer on Kinkazan also occurs because of the occasional high deer density on the island. Similarly, sika deer in Nara inhabit an urbanized environment (Nara Park), and their population is controlled by humans. Large decreases in population size have occurred in the past. In addition, these deer exhibit dwarfishness in body size and delay of sexual maturity (Torii and Tatsuzawa, Reference Torii, Tatsuzawa, MacCullough, Kaji and Takatsuki2009). These phenomena are thought to be a result of selection when the deer in Nara Park are at high density. These unusual environmental conditions and/or population bottlenecks probably resulted in accelerated morphological evolution in deer lice that is not concordant with their phylogenetic history.
ACKNOWLEDGEMENTS
We thank Koichi Ikeda, Ryo Kinomoto, Toru Koizumi, Koji Mizota, Hiroshi Mizuno, Jack Mortenson, Takahiro Ohba, Satoshi D. Ohdachi, Ayuko Ohkawa, Hiroshi Takahashi, Michihisa Takatsuka, Yoshihisa Takimoto, Shirow Tatsuzawa, Chisato Terada, Tsuneaki Yabe, members of Kamikatsu Town Hunting Club, and the Foundation for the Protection of Deer in Nara Park for help in sampling. Takeo Yamauchi provided information regarding hippoboscid flies parasitizing deer. We also thank Shin-ichi Akimoto for instruction in the statistical analyses. Comments from an anonymous referee were valuable to make our arguments clearer. This study was originally conducted as A.M.'s Master's thesis, and he thanks the staff and members of Systematic Entomology, Hokkaido University for their suggestions and encouragement.
FINANCIAL SUPPORT
This study was supported by a Grand-in-Aid from the Japan Society for the Promotion of Science Nos. 18770058 and 21770083 to K.Y., 2337003701 to K.Y., as project leader S. Akimoto (Hokkaido Univ.), and National Science Foundation grants DEB-0612938 and DEB-1050706 to K.P.J. We declare that no competing interest exists.
