Ticks and Tick-borne Diseases 5 (2014) 152–160
Contents lists available at ScienceDirect
Ticks and Tick-borne Diseases journal homepage: www.elsevier.com/locate/ttbdis
Original article
Multilocus sequence typing using mitochondrial genes (mtMLST) reveals geographic population structure of Ixodes ricinus ticks Ruth E. Dinnis a,1 , Frederik Seelig a,∗,1 , Antra Bormane b , Michael Donaghy c , Stephanie A. Vollmer a , Edward J. Feil a , Klaus Kurtenbach a , Gabriele Margos a,d,e a
University of Bath, Department of Biology and Biochemistry, Claverton Down, Bath BA2 7AY, UK State Agency Infectology Center of Latvia, Linezera Strr. 3, LV-1006 Riga, Latvia c Department of Clinical Neurology, John Radcliffe Hospital, University of Oxford, Oxford OX3 9DU, UK d National Reference Centre for Borrelia at the Bavarian Authority for Public Health and Food Safety, Branch Oberschleißheim, Veterinärstraße 2, 85764 Oberschleißheim, Germany e Institute for Infectious Diseases and Zoonoses, Ludwig-Maximilians-University Munich, Germany b
a r t i c l e
i n f o
Article history: Received 12 December 2011 Received in revised form 8 August 2013 Accepted 10 October 2013 Available online 19 December 2013 Keywords: Ixodes ricinus Biogeography Population structure Multilocus sequence typing Mitochondrial DNA
a b s t r a c t The hard tick Ixodes ricinus is the principal vector of Lyme borreliosis (LB) group spirochaetes in Europe, but it also transmits a large number of other microbial pathogens that are of importance to animal and human health. Here, we characterise geographically distinct populations of this important ectoparasite based on multilocus sequence typing (MLST) of multiple mitochondrial (mt) genes (mtMLST). Internal fragments of approximately 500 bp were amplified and sequenced for 6 protein-encoding and ribosomal genes (atp6, coi, coii, coiii, cytB, and 12s). The samples analysed consisted of 506 questing nymphs collected in Britain and Latvia in 2006–2008 and in Latvia in 2002. Although little genetic structure has previously been observed in I. ricinus ticks among Europe, our data could clearly differentiate these 2 populations. Here, we argue that this novel scheme provides additional phylogenetic resolution which is important for understanding the genetic and geographic structure of I. ricinus populations. This in turn will benefit monitoring and management of tick-borne diseases. © 2013 Elsevier GmbH. All rights reserved.
Introduction The hard tick Ixodes ricinus is the principal vector for the bacteria causing Lyme borreliosis (LB) in Europe. In addition, I. ricinus transmits a large number of other microbial and viral pathogens that are of importance to animal and human health, with the most recent addition being a relapsing fever spirochaete, Borrelia miyamotoi (Platonov et al., 2011). The distribution of I. ricinus ranges from Scandinavia to western Russia and extends to the south as far as to the Atlas Mountains in northern Africa (Blaˇskovic, 1967; Yousfi-Monod and Aeschlimann, 1986). Early studies used cuticular hydrocarbon analysis to discriminate between European strains ˜ et al., 1996, 1998) and found a certain of I. ricinus (Estrada-Pena amount of population structure. However, simultaneous and subsequent studies using 18s rDNA (Mangold et al., 1998), allozymes (Delaye et al., 1997), mitochondrial (mt) genes (Casati et al., 2008; Xu et al., 2003; Noureddine et al., 2011) and microsatellites (De
∗ Corresponding author at: Zuid-Brouwerstraat 18, 2013 WT Haarlem, The Netherlands. Tel.: +31 646396971. E-mail address: [email protected] (F. Seelig). 1 These authors contributed equally to the work. 1877-959X/$ – see front matter © 2013 Elsevier GmbH. All rights reserved. http://dx.doi.org/10.1016/j.ttbdis.2013.10.001
Meeûs et al., 2002) did not find any structure within European I. ricinus populations. Delaye et al. (1997) investigated the population structure of I. ricinus from 5 different sites in Switzerland, including foci of tick-borne encephalitis (TBE), using allozymic data [␣-glycerophosphate dehydrogenase (␣-GPD) and phosphoglucomutase (PGM)]. Despite differences in pathogen prevalence, no evidence of geographic structuring of the tick populations was found. However, the study was limited to samples from a single country (Switzerland) and used only 2 polymorphic allozymes. Casati et al. (2008) studied I. ricinus samples from various European countries using mt markers (coi, coii, cytB, 12s, and the mt control region, CR), but these authors also reported a lack of geographic structure. This auspicious approach was limited by a small sample size (12 ticks from 6 countries) and an imbalanced sequence contribution of the various loci. More recently, Noureddine et al. (2011) used a combination of 2 mt and 4 nuclear gene fragments and confirmed the absence of genetic structure among 60 I. ricinus ticks sampled across Europe. This contrasted with samples from northern Africa, which formed a divergent clade. This divergence could be the result of genetic drift following geographic isolation and/or selection pressure due to different ecological conditions. However, in this study only one tick sample from each of 40 European sites
R.E. Dinnis et al. / Ticks and Tick-borne Diseases 5 (2014) 152–160
outside its local scale in France was used, thus making the single samples less representative of these populations. Kempf et al. (2011) used 11 microsatellite markers to test for host-associated genetic structure of about 600 I. ricinus ticks collected directly from host animals captured in France, Belgium, and Slovakia. They reported significant levels of genetic structuring within tick populations and suggested that host usage is not random in I. ricinus. The aim of this study was to study the genetic diversity of ticks within small-scale populations and their association with certain host types and to a lesser extent to compare I. ricinus populations across Europe. For other Ixodes species with a wide geographic distribution range, population structure has been reported. These studies have focused on 2 Ixodes species also capable of transmitting LB group spirochaetes: I. scapularis in North America and I. uriae, which lives in seabird colonies in the Atlantic. Rich et al. (1995) and Norris et al. (1996) used the 12s and 16s mt genes to study I. scapularis populations from the East coast of the United States (US). These and additional studies suggested that 2 distinct I. scapularis populations exist, one in the north-eastern region (called the American clade) and one in the south-eastern region (termed the Southern clade) of the US (Rich et al., 1995; Norris et al., 1996; Qiu et al., 2002). In I. uriae, Kempf et al. (2009) found no evidence of geographic structuring based on the mt gene cytochrome oxidase III (coiii), whereas McCoy et al. (2003) did find population differentiation using microsatellite markers. In the latter study, it was argued that the dispersal of ticks was tightly linked to the migration of their host (in this case the black-legged kittiwake, Rissa tridactyla). However, it should be noted that I. uriae is only distantly related to the I. persulcatus-I. ricinus complex and that the organisation of its mitochondrial genome differs from that of other Ixodes species (Shao et al., 2005). For bacterial pathogens, multilocus sequence typing (MLST) has been used to great effect to infer population structure (Gevers et al., 2005). While many different methods have been investigated to monitor the spread of pathogens (reviewed in Maiden, 2006), MLST has emerged as an invaluable tool to monitor the spread of virulent bacterial strains or fungal pathogens (Enright et al., 2000; Maiden et al., 1998; Meyer et al., 2009). The advantages include reproducibility of typing methods between laboratories, portability and data sharing via the internet, and more robust analyses than single loci can offer (Maiden et al., 1998; Urwin and Maiden, 2003). In addition, MLST is based on targeted PCR of several loci, which renders it applicable to field-derived samples. There are several aspects that need to be taken into consideration when developing MLST schemes. The genes should occur in single copies and should not be prone to recombination. In prokaryotic studies, slowly evolving housekeeping genes, which are under purifying selection, have traditionally been chosen (Maiden, 2006). In general, mt markers have proven useful for investigating population structure due to high mutation and low recombination rates (Li, 1997), and sequence-based markers may have certain advantages over non-sequence-based markers (Brito and Edwards, 2009). In addition, mt genes have all of the required characteristics appropriate for MLST schemes and have been used successfully as part of a scheme studying Cryptococcus neoformans and Cryptococcus gattii (Meyer et al., 2009), the causative agents of fungal meningitis (Safdieh et al., 2008). Here, we describe the development of an MLST scheme based on mt genes to study I. ricinus populations. To increase the phylogenetic resolution of mt markers, 6 genes were chosen to develop the scheme [ATPase 6 (atp6), cytochrome oxidase I (coi), II (coii), III (coiii), the small rRNA subunit (12s), and cytochrome B (cytB)]. The scheme was tested using field-collected questing I. ricinus nymphs from 2 geographically separated regions, i.e., Latvia and Britain. Using phylogenetic and goeBURST analyses,
153
Table 1 Summary of ticks collected per year and per site in Britain and Latvia. Lowercase letters in brackets correspond to symbols on the map in Fig. 1. Country
Site
Year
No. of ticks
Total no.
American Museum (a)
2006 2006 2007 2008 2007 2008 2008 2007 2008 2006 2007 2008 2009 2006 2006 2007
6 6 19 15 4 6 28 61 28 69 12 17 39 1 1 6
318
2002 2006 2007 2002 2006 2007 2002 2006 2007
4 58 18 17 20 14 14 26 17
188
Bathampton Woods (b) Eastwood (c)
Britain
Exmoor (d) Inverness (e) New Forest (f) Rainbow Woods (g) Richmond Park (h) Thurlbear Woods (i) Warleigh (j)
Babite/Tireli (k)
Latvia
Jaunciems (l)
Jurmala (m)
506
we were able to demonstrate that ticks clustered mostly according to their geographic origin. Materials and methods Tick samples used in this study Questing nymphal and adult ticks were collected by dragging a white cotton blanket (1 m2 ) attached to a wooden dowel over the vegetation. The blanket was checked for attached ticks after every 20 steps (approximately 10 m). Ticks were removed with forceps and placed in a 1.5-ml Eppendorf tube filled with 70% ethanol. Questing ticks were identified using morphological criteria using the classification keys provided by Hillyard (1996). Due to the large numbers of questing ticks collected in Britain, a subset was randomly chosen for specific identification. Because of their ecology, it can be assumed that the vast majority of all ticks collected by blanket dragging are I. ricinus, while other species with a nidicolous (nest-adapted) lifestyle, such as I. hexagonus, are very unlikely to be picked up (Milne, 1943, cited in Hillyard, 1996). Ticks were classified according to stage, i.e., larva, nymph, or adult. I. ricinus nymphs from Latvia collected in 2002 and from 2006 to 2008 (188 samples) and Britain collected from 2006 to 2008 (318 samples) were used for the comparative analysis (Table 1). In Britain, questing nymphs were collected in temperate woodland from 6 locations near Bath, England (Fig. 1A) and from one site in Richmond Park in London, from Exmoor, the New Forest, and a site near Inverness, Scotland (see Fig. 1 for an overview of all locations). Tick samples from Latvia were collected from 3 mixed coniferous/deciduous forest sites surrounding Riga, Latvia (Fig. 1B). See Etti et al. (2003) for a detailed habitat description of the 3 Latvian sites. DNA extraction from sampled ticks Total genomic DNA from questing I. ricinus ticks was prepared using the alkaline hydrolysis method described previously (Guy
154
R.E. Dinnis et al. / Ticks and Tick-borne Diseases 5 (2014) 152–160
Fig. 1. Location of tick collection sites in Europe (A). (B) shows sites near Bath, UK. (C) shows location of sites near Riga, Latvia. Lowercase letters depict sites within a country. UK: American Museum (a); Bathampton Woods (b); Eastwood (c); Exmoor (d); Inverness (e); New Forest (f); Rainbow Woods (g); Richmond Park (h); Thurlbear Woods (i); Warleigh (j). Latvia: Babite/Tireli (k); Jaunciems (l); Jurmala (m).
and Stanek, 1991). Ticks that were partially or fully engorged with blood were extracted using ammonia hydrolysis initially, followed by DNA purification with a Qiagen DNeasy Blood, and Tissue extraction kit (Qiagen, Germany). DNA purification was carried out mostly according to the manufacturer’s instructions, apart from the 2 initial steps, in which 200 l of ATL buffer was added to the tick lysate. Thirty l of proteinase K was added to this mixture and incubated at 56 ◦ C for 12–18 h. Samples were eluted in the provided elution buffer with 2 elution steps of 100 l rather than the
suggested single step of 200 l. Extracted DNA was used directly in PCR applications. Selection of mitochondrial genes and primer design No full mt genome sequence of I. ricinus has been published yet, thus, primers could not be designed directly. We therefore used published sequences from I. persulcatus (GenBank accession no. NC 004370.1) and I. hexagonus (GenBank accession no.
R.E. Dinnis et al. / Ticks and Tick-borne Diseases 5 (2014) 152–160
NC 002010.1), which belong or are rather closely related to the Ixodes persulcatus species complex, and from I. uriae (GenBank accession no. NC 006078.1). Genes from these mt genomes were used for primer design by aligning the sequences in MEGA version 4.1 and designing primers in conserved regions. These primers are shown in Supporting Information (SI) Table 2. Genes encoding tRNA were dismissed as too small since it was our preference to have 400–500 bp of sequence data available for each gene. PCR and sequencing PCR was performed using Bioline BiomixTM master mix at 1× concentration (Bioline, UK). Primers of stock solution of 10 pmol were added to PCR reactions; volumes are shown in SI Table 3. Two l of template DNA was added per PCR reaction. The amount of template DNA was increased to 5 l if PCR reactions failed and required repeating. Reaction volumes were adjusted to 25 l with sterile distilled water. PCR samples were processed according to the PCR and optimisation criteria as detailed in SI Table 3. Using a UV transilluminator, PCR products were visualised in agarose gels. A minority of PCR reactions produced multiple bands. The main band of the expected size was cut out of the agarose gel and purified using a QIAquick kit (Qiagen, Germany). Forward and reverse nucleotide sequences of PCR amplicons were sequenced by Qiagen Genomics, Germany, and Agencourt, USA, using the same primers as for PCR amplification (SI Table 2). The quality of forward and reverse sequences was compared using DNASTAR Lasergene SeqMan (version 7.1.0). Traces were checked for errors by using an automatic ambiguous base matching functions and by manually looking at poor quality bases. DNA sequences of sufficient quality were aligned using the default setting for ClustalW alignment and trimmed according to the selected allele type region. Gene statistics Both the ratios of overall non-synonymous and synonymous substitutions (dN/dS) of concatenated sequences and the mean nucleotide diversity per site () (for all samples and for British and Latvian samples independently) were determined in MEGA version 4.1 (Tamura et al., 2007). For the dN/dS ratios, we used the modified Nei-Gojobori method and the p-distance model, while coding and non-coding domains were assigned to measure nucleotide diversity. All positions containing alignment gaps and missing data were eliminated only in pairwise sequence comparisons (pairwise deletion option). MLST scheme For each gene, alleles were assigned to sequences using a nonredundant database (NRDB) (Gish, 2008). Allele numbers were assigned to unique sequences with the first discovered allele given the type number 1, the second 2, and so on. New allele types were verified by aligning the new sequence with all existing alleles of that gene in MEGA version 4.1 (Tamura et al., 2007). Differing bases were verified by examining the original traces. Consecutive allele numbers were assigned if no other allele matched the new sequence. Alleles with identical sequences were given the same allele designation. The allelic profiles were combined in the order of atp6, coi, coii, coiii, 12s, and cytB to create a 6-integer allelic profile and to determine sequence types (ST). The allelic profiles were assigned arbitrary ST numbers, and identical combinations of alleles were given identical ST numbers. New allelic profiles were given new ST designations. After allocation of sequence types according to the allelic profile, singleton STs were compared with the entire set of sequence types
155
to find the closest match. New allele types, which differed from the more common sequence types, were reanalysed and differing bases checked once again in the sequence traces to further confirm the sequence identity. Sequence analysis We used the software Phyloviz (http://www.phyloviz.net) for MLST analysis, which employs the eBURST algorithm (Francisco et al., 2009; Feil et al., 2004). Samples with the same STs form nodes, which are proportionally sized to reflect the frequency of the ST within the population, i.e., larger nodes indicate STs that are more common. Depending on the number of single (SLV), double (DLV), or triple locus variants (TLV), nodes are related to each other with links. A permutation test (Margos et al., 2008) was used to test whether genetic differences between populations were significant. Arlequin statistical analyses The program Arlequin version 3.1 (Excoffier et al., 2005) was used to investigate population differentiation and to perform an analysis of molecular variance (AMOVA). Populations were created using the grouping function according to the country of origin. PhiST values were computed at a significance level of 0.05 with 10,000 permutations. Values of PhiST can range from zero to one. Values of zero for PhiST indicate that the populations are completely homogeneous, values of one indicate disparate populations, values greater than 0.25 indicate a large amount of genetic diversity (Freeland, 2005). Phylogenetic analyses Alignments of sequences were tested for suitable models in FindModel (Tao et al., 2009). FindModel incorporates several processes to determine the best model for the submitted data. Weighbour trees (Bruno et al., 2000) based on Jukes-Cantor distances were used as starting trees. PAML (Yang, 1997, 2007) was used to calculate likelihoods. Akaike information criterion (AIC) scores were calculated using MODELTEST (Posada and Crandall, 1998). Maximum likelihood trees were constructed using PhyML 3.0 (Guindon and Gascuel, 2003) hosted at the ATGC Montpellier bioinformatics platform. The substitution model was determined in MODELTEST for each data set. Starting trees were set as BIONJ with tree improvement setting of Sub-tree Pruning and Regrafting (SPR) and Nearest Neighbour Interchange (NNI) and approximate Likelihood Ratio Test (aLRT) SH (Shimodaira-Hasegawa)-like branch support parameters. All other parameters were set as default values. Results Mitochondrial gene statistics Internal fragments of approximately 450–500 bp of 6 loci (atp6, coi, coii, coiii, 12s, and cytB) were sequenced and used to develop the mtMLST scheme. The characteristics of these loci are shown in Table 2. These data indicate that the nucleotide diversity per site () was similar in the 2 populations. Mean values of concatenated samples from Britain and Latvia were 0.0042 and 0.0059, respectively. Observed dN/dS values of
Contents lists available at ScienceDirect
Ticks and Tick-borne Diseases journal homepage: www.elsevier.com/locate/ttbdis
Original article
Multilocus sequence typing using mitochondrial genes (mtMLST) reveals geographic population structure of Ixodes ricinus ticks Ruth E. Dinnis a,1 , Frederik Seelig a,∗,1 , Antra Bormane b , Michael Donaghy c , Stephanie A. Vollmer a , Edward J. Feil a , Klaus Kurtenbach a , Gabriele Margos a,d,e a
University of Bath, Department of Biology and Biochemistry, Claverton Down, Bath BA2 7AY, UK State Agency Infectology Center of Latvia, Linezera Strr. 3, LV-1006 Riga, Latvia c Department of Clinical Neurology, John Radcliffe Hospital, University of Oxford, Oxford OX3 9DU, UK d National Reference Centre for Borrelia at the Bavarian Authority for Public Health and Food Safety, Branch Oberschleißheim, Veterinärstraße 2, 85764 Oberschleißheim, Germany e Institute for Infectious Diseases and Zoonoses, Ludwig-Maximilians-University Munich, Germany b
a r t i c l e
i n f o
Article history: Received 12 December 2011 Received in revised form 8 August 2013 Accepted 10 October 2013 Available online 19 December 2013 Keywords: Ixodes ricinus Biogeography Population structure Multilocus sequence typing Mitochondrial DNA
a b s t r a c t The hard tick Ixodes ricinus is the principal vector of Lyme borreliosis (LB) group spirochaetes in Europe, but it also transmits a large number of other microbial pathogens that are of importance to animal and human health. Here, we characterise geographically distinct populations of this important ectoparasite based on multilocus sequence typing (MLST) of multiple mitochondrial (mt) genes (mtMLST). Internal fragments of approximately 500 bp were amplified and sequenced for 6 protein-encoding and ribosomal genes (atp6, coi, coii, coiii, cytB, and 12s). The samples analysed consisted of 506 questing nymphs collected in Britain and Latvia in 2006–2008 and in Latvia in 2002. Although little genetic structure has previously been observed in I. ricinus ticks among Europe, our data could clearly differentiate these 2 populations. Here, we argue that this novel scheme provides additional phylogenetic resolution which is important for understanding the genetic and geographic structure of I. ricinus populations. This in turn will benefit monitoring and management of tick-borne diseases. © 2013 Elsevier GmbH. All rights reserved.
Introduction The hard tick Ixodes ricinus is the principal vector for the bacteria causing Lyme borreliosis (LB) in Europe. In addition, I. ricinus transmits a large number of other microbial and viral pathogens that are of importance to animal and human health, with the most recent addition being a relapsing fever spirochaete, Borrelia miyamotoi (Platonov et al., 2011). The distribution of I. ricinus ranges from Scandinavia to western Russia and extends to the south as far as to the Atlas Mountains in northern Africa (Blaˇskovic, 1967; Yousfi-Monod and Aeschlimann, 1986). Early studies used cuticular hydrocarbon analysis to discriminate between European strains ˜ et al., 1996, 1998) and found a certain of I. ricinus (Estrada-Pena amount of population structure. However, simultaneous and subsequent studies using 18s rDNA (Mangold et al., 1998), allozymes (Delaye et al., 1997), mitochondrial (mt) genes (Casati et al., 2008; Xu et al., 2003; Noureddine et al., 2011) and microsatellites (De
∗ Corresponding author at: Zuid-Brouwerstraat 18, 2013 WT Haarlem, The Netherlands. Tel.: +31 646396971. E-mail address: [email protected] (F. Seelig). 1 These authors contributed equally to the work. 1877-959X/$ – see front matter © 2013 Elsevier GmbH. All rights reserved. http://dx.doi.org/10.1016/j.ttbdis.2013.10.001
Meeûs et al., 2002) did not find any structure within European I. ricinus populations. Delaye et al. (1997) investigated the population structure of I. ricinus from 5 different sites in Switzerland, including foci of tick-borne encephalitis (TBE), using allozymic data [␣-glycerophosphate dehydrogenase (␣-GPD) and phosphoglucomutase (PGM)]. Despite differences in pathogen prevalence, no evidence of geographic structuring of the tick populations was found. However, the study was limited to samples from a single country (Switzerland) and used only 2 polymorphic allozymes. Casati et al. (2008) studied I. ricinus samples from various European countries using mt markers (coi, coii, cytB, 12s, and the mt control region, CR), but these authors also reported a lack of geographic structure. This auspicious approach was limited by a small sample size (12 ticks from 6 countries) and an imbalanced sequence contribution of the various loci. More recently, Noureddine et al. (2011) used a combination of 2 mt and 4 nuclear gene fragments and confirmed the absence of genetic structure among 60 I. ricinus ticks sampled across Europe. This contrasted with samples from northern Africa, which formed a divergent clade. This divergence could be the result of genetic drift following geographic isolation and/or selection pressure due to different ecological conditions. However, in this study only one tick sample from each of 40 European sites
R.E. Dinnis et al. / Ticks and Tick-borne Diseases 5 (2014) 152–160
outside its local scale in France was used, thus making the single samples less representative of these populations. Kempf et al. (2011) used 11 microsatellite markers to test for host-associated genetic structure of about 600 I. ricinus ticks collected directly from host animals captured in France, Belgium, and Slovakia. They reported significant levels of genetic structuring within tick populations and suggested that host usage is not random in I. ricinus. The aim of this study was to study the genetic diversity of ticks within small-scale populations and their association with certain host types and to a lesser extent to compare I. ricinus populations across Europe. For other Ixodes species with a wide geographic distribution range, population structure has been reported. These studies have focused on 2 Ixodes species also capable of transmitting LB group spirochaetes: I. scapularis in North America and I. uriae, which lives in seabird colonies in the Atlantic. Rich et al. (1995) and Norris et al. (1996) used the 12s and 16s mt genes to study I. scapularis populations from the East coast of the United States (US). These and additional studies suggested that 2 distinct I. scapularis populations exist, one in the north-eastern region (called the American clade) and one in the south-eastern region (termed the Southern clade) of the US (Rich et al., 1995; Norris et al., 1996; Qiu et al., 2002). In I. uriae, Kempf et al. (2009) found no evidence of geographic structuring based on the mt gene cytochrome oxidase III (coiii), whereas McCoy et al. (2003) did find population differentiation using microsatellite markers. In the latter study, it was argued that the dispersal of ticks was tightly linked to the migration of their host (in this case the black-legged kittiwake, Rissa tridactyla). However, it should be noted that I. uriae is only distantly related to the I. persulcatus-I. ricinus complex and that the organisation of its mitochondrial genome differs from that of other Ixodes species (Shao et al., 2005). For bacterial pathogens, multilocus sequence typing (MLST) has been used to great effect to infer population structure (Gevers et al., 2005). While many different methods have been investigated to monitor the spread of pathogens (reviewed in Maiden, 2006), MLST has emerged as an invaluable tool to monitor the spread of virulent bacterial strains or fungal pathogens (Enright et al., 2000; Maiden et al., 1998; Meyer et al., 2009). The advantages include reproducibility of typing methods between laboratories, portability and data sharing via the internet, and more robust analyses than single loci can offer (Maiden et al., 1998; Urwin and Maiden, 2003). In addition, MLST is based on targeted PCR of several loci, which renders it applicable to field-derived samples. There are several aspects that need to be taken into consideration when developing MLST schemes. The genes should occur in single copies and should not be prone to recombination. In prokaryotic studies, slowly evolving housekeeping genes, which are under purifying selection, have traditionally been chosen (Maiden, 2006). In general, mt markers have proven useful for investigating population structure due to high mutation and low recombination rates (Li, 1997), and sequence-based markers may have certain advantages over non-sequence-based markers (Brito and Edwards, 2009). In addition, mt genes have all of the required characteristics appropriate for MLST schemes and have been used successfully as part of a scheme studying Cryptococcus neoformans and Cryptococcus gattii (Meyer et al., 2009), the causative agents of fungal meningitis (Safdieh et al., 2008). Here, we describe the development of an MLST scheme based on mt genes to study I. ricinus populations. To increase the phylogenetic resolution of mt markers, 6 genes were chosen to develop the scheme [ATPase 6 (atp6), cytochrome oxidase I (coi), II (coii), III (coiii), the small rRNA subunit (12s), and cytochrome B (cytB)]. The scheme was tested using field-collected questing I. ricinus nymphs from 2 geographically separated regions, i.e., Latvia and Britain. Using phylogenetic and goeBURST analyses,
153
Table 1 Summary of ticks collected per year and per site in Britain and Latvia. Lowercase letters in brackets correspond to symbols on the map in Fig. 1. Country
Site
Year
No. of ticks
Total no.
American Museum (a)
2006 2006 2007 2008 2007 2008 2008 2007 2008 2006 2007 2008 2009 2006 2006 2007
6 6 19 15 4 6 28 61 28 69 12 17 39 1 1 6
318
2002 2006 2007 2002 2006 2007 2002 2006 2007
4 58 18 17 20 14 14 26 17
188
Bathampton Woods (b) Eastwood (c)
Britain
Exmoor (d) Inverness (e) New Forest (f) Rainbow Woods (g) Richmond Park (h) Thurlbear Woods (i) Warleigh (j)
Babite/Tireli (k)
Latvia
Jaunciems (l)
Jurmala (m)
506
we were able to demonstrate that ticks clustered mostly according to their geographic origin. Materials and methods Tick samples used in this study Questing nymphal and adult ticks were collected by dragging a white cotton blanket (1 m2 ) attached to a wooden dowel over the vegetation. The blanket was checked for attached ticks after every 20 steps (approximately 10 m). Ticks were removed with forceps and placed in a 1.5-ml Eppendorf tube filled with 70% ethanol. Questing ticks were identified using morphological criteria using the classification keys provided by Hillyard (1996). Due to the large numbers of questing ticks collected in Britain, a subset was randomly chosen for specific identification. Because of their ecology, it can be assumed that the vast majority of all ticks collected by blanket dragging are I. ricinus, while other species with a nidicolous (nest-adapted) lifestyle, such as I. hexagonus, are very unlikely to be picked up (Milne, 1943, cited in Hillyard, 1996). Ticks were classified according to stage, i.e., larva, nymph, or adult. I. ricinus nymphs from Latvia collected in 2002 and from 2006 to 2008 (188 samples) and Britain collected from 2006 to 2008 (318 samples) were used for the comparative analysis (Table 1). In Britain, questing nymphs were collected in temperate woodland from 6 locations near Bath, England (Fig. 1A) and from one site in Richmond Park in London, from Exmoor, the New Forest, and a site near Inverness, Scotland (see Fig. 1 for an overview of all locations). Tick samples from Latvia were collected from 3 mixed coniferous/deciduous forest sites surrounding Riga, Latvia (Fig. 1B). See Etti et al. (2003) for a detailed habitat description of the 3 Latvian sites. DNA extraction from sampled ticks Total genomic DNA from questing I. ricinus ticks was prepared using the alkaline hydrolysis method described previously (Guy
154
R.E. Dinnis et al. / Ticks and Tick-borne Diseases 5 (2014) 152–160
Fig. 1. Location of tick collection sites in Europe (A). (B) shows sites near Bath, UK. (C) shows location of sites near Riga, Latvia. Lowercase letters depict sites within a country. UK: American Museum (a); Bathampton Woods (b); Eastwood (c); Exmoor (d); Inverness (e); New Forest (f); Rainbow Woods (g); Richmond Park (h); Thurlbear Woods (i); Warleigh (j). Latvia: Babite/Tireli (k); Jaunciems (l); Jurmala (m).
and Stanek, 1991). Ticks that were partially or fully engorged with blood were extracted using ammonia hydrolysis initially, followed by DNA purification with a Qiagen DNeasy Blood, and Tissue extraction kit (Qiagen, Germany). DNA purification was carried out mostly according to the manufacturer’s instructions, apart from the 2 initial steps, in which 200 l of ATL buffer was added to the tick lysate. Thirty l of proteinase K was added to this mixture and incubated at 56 ◦ C for 12–18 h. Samples were eluted in the provided elution buffer with 2 elution steps of 100 l rather than the
suggested single step of 200 l. Extracted DNA was used directly in PCR applications. Selection of mitochondrial genes and primer design No full mt genome sequence of I. ricinus has been published yet, thus, primers could not be designed directly. We therefore used published sequences from I. persulcatus (GenBank accession no. NC 004370.1) and I. hexagonus (GenBank accession no.
R.E. Dinnis et al. / Ticks and Tick-borne Diseases 5 (2014) 152–160
NC 002010.1), which belong or are rather closely related to the Ixodes persulcatus species complex, and from I. uriae (GenBank accession no. NC 006078.1). Genes from these mt genomes were used for primer design by aligning the sequences in MEGA version 4.1 and designing primers in conserved regions. These primers are shown in Supporting Information (SI) Table 2. Genes encoding tRNA were dismissed as too small since it was our preference to have 400–500 bp of sequence data available for each gene. PCR and sequencing PCR was performed using Bioline BiomixTM master mix at 1× concentration (Bioline, UK). Primers of stock solution of 10 pmol were added to PCR reactions; volumes are shown in SI Table 3. Two l of template DNA was added per PCR reaction. The amount of template DNA was increased to 5 l if PCR reactions failed and required repeating. Reaction volumes were adjusted to 25 l with sterile distilled water. PCR samples were processed according to the PCR and optimisation criteria as detailed in SI Table 3. Using a UV transilluminator, PCR products were visualised in agarose gels. A minority of PCR reactions produced multiple bands. The main band of the expected size was cut out of the agarose gel and purified using a QIAquick kit (Qiagen, Germany). Forward and reverse nucleotide sequences of PCR amplicons were sequenced by Qiagen Genomics, Germany, and Agencourt, USA, using the same primers as for PCR amplification (SI Table 2). The quality of forward and reverse sequences was compared using DNASTAR Lasergene SeqMan (version 7.1.0). Traces were checked for errors by using an automatic ambiguous base matching functions and by manually looking at poor quality bases. DNA sequences of sufficient quality were aligned using the default setting for ClustalW alignment and trimmed according to the selected allele type region. Gene statistics Both the ratios of overall non-synonymous and synonymous substitutions (dN/dS) of concatenated sequences and the mean nucleotide diversity per site () (for all samples and for British and Latvian samples independently) were determined in MEGA version 4.1 (Tamura et al., 2007). For the dN/dS ratios, we used the modified Nei-Gojobori method and the p-distance model, while coding and non-coding domains were assigned to measure nucleotide diversity. All positions containing alignment gaps and missing data were eliminated only in pairwise sequence comparisons (pairwise deletion option). MLST scheme For each gene, alleles were assigned to sequences using a nonredundant database (NRDB) (Gish, 2008). Allele numbers were assigned to unique sequences with the first discovered allele given the type number 1, the second 2, and so on. New allele types were verified by aligning the new sequence with all existing alleles of that gene in MEGA version 4.1 (Tamura et al., 2007). Differing bases were verified by examining the original traces. Consecutive allele numbers were assigned if no other allele matched the new sequence. Alleles with identical sequences were given the same allele designation. The allelic profiles were combined in the order of atp6, coi, coii, coiii, 12s, and cytB to create a 6-integer allelic profile and to determine sequence types (ST). The allelic profiles were assigned arbitrary ST numbers, and identical combinations of alleles were given identical ST numbers. New allelic profiles were given new ST designations. After allocation of sequence types according to the allelic profile, singleton STs were compared with the entire set of sequence types
155
to find the closest match. New allele types, which differed from the more common sequence types, were reanalysed and differing bases checked once again in the sequence traces to further confirm the sequence identity. Sequence analysis We used the software Phyloviz (http://www.phyloviz.net) for MLST analysis, which employs the eBURST algorithm (Francisco et al., 2009; Feil et al., 2004). Samples with the same STs form nodes, which are proportionally sized to reflect the frequency of the ST within the population, i.e., larger nodes indicate STs that are more common. Depending on the number of single (SLV), double (DLV), or triple locus variants (TLV), nodes are related to each other with links. A permutation test (Margos et al., 2008) was used to test whether genetic differences between populations were significant. Arlequin statistical analyses The program Arlequin version 3.1 (Excoffier et al., 2005) was used to investigate population differentiation and to perform an analysis of molecular variance (AMOVA). Populations were created using the grouping function according to the country of origin. PhiST values were computed at a significance level of 0.05 with 10,000 permutations. Values of PhiST can range from zero to one. Values of zero for PhiST indicate that the populations are completely homogeneous, values of one indicate disparate populations, values greater than 0.25 indicate a large amount of genetic diversity (Freeland, 2005). Phylogenetic analyses Alignments of sequences were tested for suitable models in FindModel (Tao et al., 2009). FindModel incorporates several processes to determine the best model for the submitted data. Weighbour trees (Bruno et al., 2000) based on Jukes-Cantor distances were used as starting trees. PAML (Yang, 1997, 2007) was used to calculate likelihoods. Akaike information criterion (AIC) scores were calculated using MODELTEST (Posada and Crandall, 1998). Maximum likelihood trees were constructed using PhyML 3.0 (Guindon and Gascuel, 2003) hosted at the ATGC Montpellier bioinformatics platform. The substitution model was determined in MODELTEST for each data set. Starting trees were set as BIONJ with tree improvement setting of Sub-tree Pruning and Regrafting (SPR) and Nearest Neighbour Interchange (NNI) and approximate Likelihood Ratio Test (aLRT) SH (Shimodaira-Hasegawa)-like branch support parameters. All other parameters were set as default values. Results Mitochondrial gene statistics Internal fragments of approximately 450–500 bp of 6 loci (atp6, coi, coii, coiii, 12s, and cytB) were sequenced and used to develop the mtMLST scheme. The characteristics of these loci are shown in Table 2. These data indicate that the nucleotide diversity per site () was similar in the 2 populations. Mean values of concatenated samples from Britain and Latvia were 0.0042 and 0.0059, respectively. Observed dN/dS values of
