Mitochondrial phylogeography and population structure of the cattle tick Rhipicephalus appendiculatus in the African Great Lakes region

Mitochondrial phylogeography and population structure of the cattle tick Rhipicephalus... Background: The ixodid tick Rhipicephalus appendiculatus is the main vector of Theileria parva, wich causes the highly fatal cattle disease East Coast fever (ECF) in sub-Saharan Africa. Rhipicephalus appendiculatus populations differ in their ecology, diapause behaviour and vector competence. Thus, their expansion in new areas may change the genetic structure and consequently affect the vector-pathogen system and disease outcomes. In this study we investigated the genetic distribution of R. appendiculatus across agro-ecological zones (AEZs) in the African Great Lakes region to better understand the epidemiology of ECF and elucidate R. appendiculatus evolutionary history and biogeographical colonization in Africa. Methods: Sequencing was performed on two mitochondrial genes (cox1 and 12S rRNA) of 218 ticks collected from cattle across six AEZs along an altitudinal gradient in the Democratic Republic of Congo, Rwanda, Burundi and Tanzania. Phylogenetic relationships between tick populations were determined and evolutionary population dynamics models were assessed by mismach distribution. Results: Population genetic analysis yielded 22 cox1 and 9 12S haplotypes in a total of 209 and 126 nucleotide sequences, respectively. Phylogenetic algorithms grouped these haplotypes for both genes into two major clades (lineages A and B). We observed significant genetic variation segregating the two lineages and low structure among populations with high degree of migration. The observed high gene flow indicates population admixture between AEZs. However, reduced number of migrants was observed between lowlands and highlands. Mismatch analysis detected a signature of rapid demographic and range expansion of lineage A. The star-like pattern of isolated and published haplotypes indicates that the two lineages evolve independently and have been subjected to expansion across Africa. (Continued on next page) * Correspondence: gastonamzati@gmail.com Unit of Integrated Veterinary Research, Department of Veterinary Medicine, Faculty of Sciences, Namur Research Institute for Life Sciences (NARILIS), University of Namur (UNamur), Rue de Bruxelles 61, 5000 Namur, Belgium Research Unit of Veterinary Epidemiology and Biostatistics, Faculty of Agricultural and Environmental Sciences, Université Evangélique en Afrique, P.O. Box 3323, Bukavu, Democratic Republic of the Congo Full list of author information is available at the end of the article © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Amzati et al. Parasites & Vectors (2018) 11:329 Page 2 of 19 (Continued from previous page) Conclusions: Two sympatric R. appendiculatus lineages occur in the Great Lakes region. Lineage A, the most diverse and ubiquitous, has experienced rapid population growth and range expansion in all AEZs probably through cattle movement, whereas lineage B, the less abundant, has probably established a founder population from recent colonization events and its occurrence decreases with altitude. These two lineages are sympatric in central and eastern Africa and allopatric in southern Africa. The observed colonization pattern may strongly affect the transmission system and may explain ECF endemic instability in the tick distribution fringes. Keywords: East Coast fever, Theileria parva, Rhipicephalus appendiculatus, Phylogenetic, cox1, 12S rRNA, Ticks, Evolutionary history, Agro-ecological zones, Population genetics Background The adaptive mechanism of R. appendiculatus to chan- The ixodid brown ear tick Rhipicephalus appendiculatus ged environment suggests genetic divergence between is the main vector of the protozoan pathogen Theileria geographical stocks and phenotypic variations, including parva, the causative agent of a fatal lymphoproliferative diapause behaviour and vector competence to transmit T. cattle disease known as East Coast fever (ECF). East Coast parva [23, 24]. The degree to which different R. appendi- fever is a highly pathogenic and the most economically culatus stocks acquire and transmit T. parva is partially important tick-borne disease of cattle in 12 sub-Saharan genetically dependent because of the heritability of their African countries, including Burundi, Democratic Repub- susceptibility to infection [25]. Furthermore, there is a lic of Congo and Rwanda [1]. Rhipicephalus appendicula- similar extensive genetic variation among T. parva strains tus is the most abundant tick in the Great Lakes region of in the field associated with different clinical features and Central Africa, where its burden and distribution vary sig- disease outcomes, and variable cross-immunity levels. nificantly among agro-ecological zones (AEZs) [2–4]. The Studies suggest that there are also specific interactions geographical dispersal pattern and population dynamics of between R. appendiculatus and T. parva stocks in the this tick are mainly driven by climatic conditions, vegeta- transmission dynamic system, significantly affecting the tion, host availability and mobility, grazing system and epidemiology of ECF [26]. Thus, phylogenetic and management practices [5, 6]. ecological analyses should provide useful information to The Great Lakes region of Central Africa is charac- control ECF by determining: (i) the genetic structure of R. terised by cattle movement within and between countries appendiculatus populations; (ii) its dispersal pattern; and for trade, breeding and pasture [7, 8]. During the (iii) its demographic history [27, 28]. pre-colonial period, immigrants originating from Rwanda The genetic diversity of R. appendiculatus has been and Burundi settled with their cattle in Congo in search of studied in different African countries using various mo- grazing lands. In addition, political unrest during the past lecular tools, such as mitochondrial DNA [14, 29] and three decades and rapidly increasing demand for animal micro- and minisatellite DNA markers [30]. Two distinct products increased significantly the importation of live an- genetic groups have been described, namely the eastern imals [9, 10]. This cross-border movement of cattle across and the southern African lineages [29]. More recently, geographical areas, together with bioclimatic conditions Kanduma et al. [31] found that the two genetic groups suitable for R. appendiculatus, could play a significant role are present in Kenya. These evidences show that R. in spreading ticks and pathogens [11–14]. Therefore, the appendiculatus genetic groups may have a wide geo- spread and establishment of ticks from one geographical graphical range, with different ecological preferences region to another might be setting up a complexity in the and phenology in sub-Saharan Africa [32–34], due to epidemiological status and control of the disease they differences in body size [35] and diapause induction and transmit [15–17]. Thus, predicting vector-borne pathogen intensity [36, 37]. dynamics and emergence relies on better understanding Major gaps in knowledge still exist concerning the of mechanisms underlying the genetic structure of their agro-ecological colonization and establishment of the R. vectors [18, 19]. appendiculatus lineages in the Great Lakes region of Ecological establishment and population genetic struc- Central Africa, where cattle mobility seems to be the main ture of ticks can be affected by founder events and gene factor of tick dispersal and epidemic instability of ECF [2, flow, largely due to their dispersal across geographical 38]. Thus, further studies on the population structure and zones through host migration [18, 20]. Arthropod vectors phylogeography of R. appendiculatus, including ticks from then respond differently to evolutionary forces such as mi- distinct agro-ecological conditions of DR Congo, Burundi gration, mutation, selection and genetic drift (promoted and Rwanda, are important to shed light on the intra and by bottlenecks) in their new environment [21, 22]. inter population variation, the dispersal pattern and the Amzati et al. Parasites & Vectors (2018) 11:329 Page 3 of 19 historical dynamics of the characterised lineages in various close to the equator, the wide range of altitudes (700 to ecological situations of Africa. The objective of this study 3000 m above sea level) mitigates significantly the was to assess the evolutionary relationship between R. tropical climate attributes and therefore the AEZs limits. appendiculatus populations and to investigate the impact Generally, the rainfall period is long and bimodal (with of geographical locations on its genetic structure, to better variations in the duration between AEZs). An early rainy understand the epidemiology of ECF in the Great Lakes season starts approximately in September and ends in region. To achieve this objective, we sequenced and ana- December, while the late rainy season lasts from February lysed fragments of the cytochrome c oxidase subunit 1 to May, followed by a dry season from June to August. (cox1) and the 12S ribosomal RNA (12S rRNA) gene loci. There is an intervening short dry period between the two The genetic structure of R. appendiculatus provides clues rainy seasons of about 15 days in January and/or February to a better understanding of the epidemiology of ECF and depending on the AEZs. Rainfall and temperature are insight into the development of targeted selective manage- strongly influenced by altitude ranges. Rainfall increases ment and control strategies. while temperature decreases with increasing elevation. The average temperature ranges between 12–25 °C and Methods the annual average rainfall from 800–1100 mm in the Study area Ruzizi Valley to 1300–2000 mm in the highlands. Ticks were collected from cattle in six agro-ecological In eastern DRC, ticks were collected along an altitudinal zones (AEZs) of the Central African Great Lakes region transect in South-Kivu province. Based on elevation and (Democratic Republic of Congo, Rwanda and Burundi) geographical position within the transect, three major AEZs (Fig. 1). The three countries share the Ruzizi Valley, were defined from lowlands to highlands, namely, DRC which consists of lowlands along the Ruzizi River flow- AEZ1 (lowlands) in the Ruzizi Valley, DRC AEZ2 (mid- ing between Lake Kivu and Lake Tanganyika. A sum- lands) in the administrative district of Walungu and DRC mary description of the characteristics of the six AEZs is AEZ3 (highlands) in the district of Kabare and the high- presented in Table 1. Briefly, although the study area is lands part of Walungu (Mulumemunene). The lowlands Fig. 1 Sampling sites of Rhipicephalus appendiculatus ticks in DRC, Rwanda and Burundi. a Map of Africa showing the study area (in grey) and other countries where the tick was previously sequenced (indicated by their names). b Sampling localities of R. appendiculatus and their altitudes (squares: AEZ1 altitude < 1200 m; circles: AEZ2 altitude 1200–1600 m; and triangles: AEZ3 altitude > 1600 m). The sites represented by empty circles and triangles show sampling locations described by Mtambo et al. [29] Amzati et al. Parasites & Vectors (2018) 11:329 Page 4 of 19 Table 1 Geographical and climatic attributes of the six agro-ecological zones (AEZ) Country Agro-ecological zone (AEZ) Altitude (m) Temperature (°C) Rainfall (mm/year) Rainy season Sample size (no. of ticks) DRC AEZ1 780–1100 23–25 800–1000 October-April 46 AEZ2 1200–1600 17–21 1000–1500 September-May 54 AEZ3 1600–2800 12–19 1350–2000 September-May 46 Burundi AEZ1 774–1100 23–24 800–1100 November-May 26 AEZ3 1700–2800 14–15 1300–2000 September-May 17 Rwanda AEZ2 1200–1500 21–24 800–950 November-May 20 Notes: DRC (AEZ1: Lowlands, AEZ2: Midlands, AEZ3: Highlands); Burundi (AEZ1: Lowlands, AEZ3: Highlands); Rwanda (AEZ2: eastern low plateau which is the lowlands of Rwanda as described by Bazarusanga et al. [3]) Sequences previously described in Mtambo at al. [29] are not included in the 20 samples from Rwanda and thus were not used in the population genetic analysis presented in Tables 2, 3 and 4 region (DRC AEZ1) is the warmest AEZ, characterised by a obtained from the Sokoine University of Agriculture in tropical dry climate (semi-arid), with a cool dry season of Tanzania. 4–5 months. Cattle are raised generally in an open grazing system in communal pastures of savannah grassland. In DNA extraction, PCR amplification and sequencing contrast, the “upland Kivu” (DRC AEZ2 and AEZ3) has a Total genomic DNA was isolated from whole individual montane humid tropical climate and is much cooler. The ticks using the DNeasy® Blood and Tissue Kit (Qiagen warm dry season lasts for 3–4 months, with occasional and GmbH, Hilden, Germany) according to the standard poorly distributed rainfall. manufacturer’s protocol, except that an additional incu- Ticks from Burundi were collected from two AEZs. bation of 10 min at 56 °C was added after mixing the The sampled administrative districts included Rugombo sample with 200 μl buffer AL to ensure optimal cell lysis. and Gihanga in the Imbo lowlands (Burundi AEZ1) and Prior to extract DNA, ticks were washed twice in double Muramvya, Mwaro and Mugamba districts of the distilled water and left to dry for 10 min at room Congo-Nile highlands (Burundi AEZ3). The Imbo region temperature and homogenized by grinding. extends along the Ruzizi River and the North of Lake Given the suitability of mitochondrial genes to dis- Tanganyika. The vegetation is composed of savannah criminate intraspecific variation in ticks [29, 31], we and small patches of forest. amplified the cox1and 12S rRNA gene loci to assess the Ticks from Rwanda were obtained from the eastern low genetic diversity and phylogenetic relationships of R. plateau in Nyagatare, Gatsibo and Kayonza districts appendiculatus populations. A 710 bp fragment of cox1 (Rwanda AEZ2). This region includes the savannah of gene locus was amplified with the forward primer eastern province of Rwanda. The eastern semi-arid plateau LCO1490 (5'-GGT CAA CAA ATC ATA AAG ATA zone is the most important pastoral region in Rwanda, TTG G-3') and the reverse primer HC02198 (5'-TAA holding 40% of the national herd raised in a communal ACT TCA GGG TGA CCA AAA AAT CA-3') as previ- grazing system. The vegetation type is largely savannah ously described by Folmer et al. [40]; for the 12S rRNA and some river bank woods. gene locus we used the primers SR-J-1499 (5'-TAC TAT GTT ACG ACT TAT-3') and SR-N-14594 (5'-AAA CTA Tick sampling and morphological identification GGA TTA GAT ACC C-3') with a fragment size of 380 Ticks were collected from cattle during a cross-sectional bp, described in Simon et al. [41]. PCR amplifications of survey conducted from February to April 2015 (late both cox1 and 12S rRNA gene fragments were per- rainy season) in the six AEZs. In each AEZ, 8 to 12 cat- formed using 50 ng of genomic DNA, 25 μl of 2× Accu- tle herds were selected randomly using a random num- Power® PCR PreMix (Bioneer PCR-PreMix, Seoul, South ber generator in Microsoft Excel. From ticks collected in Korea), 0.2 μM each of forward and reverse primers, and each herd representing a location or village, 4 to 5 ticks nuclease free water added up to a final reaction volume were selected using the same random process for further of 50 μl. The thermal cycling program for cox1 consisted analysis. The number of ticks sampled in each popula- of an initial denaturation at 95 °C for 3 min followed by tion is shown in Table 1. All ticks were directly 35 cycles of denaturation at 94 °C for 1 min, annealing immersed in 70% ethanol for preservation and morpho- at 40 °C for 1 min and extension at 72 °C for 1 min. The logically identified using the identification key of Walker final extension was carried out for 10 min at 72 °C. PCR et al. [39]. Morphological identification was confirmed at parameters of 12S rRNA gene fragment were as de- the tick unit at the International Livestock Research In- scribed for cox1, except that the annealing temperature stitute (ILRI, Kenya). Ten additional specimens originat- was 52 °C and the extension time was 90 s. PCR prod- ing from Simanjiro district in northern Tanzania were ucts were analysed by electrophoresis on a 1.8% agarose Amzati et al. Parasites & Vectors (2018) 11:329 Page 5 of 19 gel. Amplicons were purified using the QIAquick® PCR migrants. These analyses were performed for combined Purification Kit (Qiagen GmbH, Hilden, Germany) fol- data (of the main haplogroups identified) to understand lowing the manufacturer’s instructions. Both strands the effect of coexistence of the haplogroups in the gen- were sequenced directly with the same primers used for etic structure and diversity. PCR, on an ABI 3730 capillary sequencer (Applied Bio- systems, California, USA). Demographic and spatial expansion history analyses The historical population dynamics and structure was Sequences editing, blasting and multiple alignments inferred from mismatch distribution of cox1 haplotypes Forward and reverse chromatograms for each individual implemented in Arlequin [45], which compared the ob- tick were visually checked. Sequences were manually served distribution of pairwise nucleotide differences be- edited and assembled using CLC Main Workbench soft- tween haplotypes and that expected under a population ware v7.8.1 (CLC Bio, Aarhus, Denmark). Multiple se- expansion model for each population, haplogroup and quences were aligned with CrustalW algorithm using overall data. Multimodal mismatch pattern is assumed default settings in the same software. The sequences to define a population of demographic equilibrium or were then trimmed to exclude poor quality bases and constant size, whereas sudden or expanding population obtain uniform sizes. The final sequence sizes were 586 is characterised by unimodal distribution. The sum of bp and 346 bp for cox1 and 12S rRNA genes, respect- square deviation (SSD) were estimated to determine if ively. Aligned and trimmed sequences were subjected to the observed mismatch deviated significantly from a a BLAST search against all NCBI nucleotide databases, population expansion model, and the Harpending’s rag- with default settings to confirm their species identity. A gedness index (RI) were used to evaluate the smoothness total of 219 sequences for cox1 and 126 for 12S were ob- of mismatch distribution [46]. In addition, to detect de- tained after quality check processing. To check for amp- viation from neutrality expectations or mutation-drift lification of nuclear pseudogenes and confirm protein equilibrium we performed analysis of Fu’s Fs [47] and coding, all cox1 nucleotide sequences were translated Tajima’s D [48] statistics in Arlequin and DnaSP. Taji- into their amino acid sequences to examine the presence ma’s D test is based on the difference between the num- of ambiguous stop codons for correct coding of inverte- ber of polymorphic sites and the mean number of brate mitochondrial DNA. nucleotide pairwise differences, while Fu’s Fs test is based on haplotype frequencies. The significance of all Genetic diversity and population genetic structure these statistics was tested by bootstrap resampling of Multiple sequences extracted from the alignments were 1000 coalescent simulations. Significant negative values used to construct haplotypes in DnaSP software v5.10.01 of neutrality statistics should indicate a signature of his- [42]. Genetic variation within populations was estimated torical event of population expansion, whereas signifi- on cox1 gene sequences. Population genetic indices rep- cantly positive values indicate events such as recent resented by number of haplotypes (h), number of segre- population bottleneck, population subdivision or pres- gating sites (S), haplotype diversity (Hd), mean number ence of some recent immigrants in a population. Values of pairwise nucleotide differences within population (K) near zero and not significant, indicate that population and nucleotide diversity (π) were calculated for each size is constant. AEZs (named populations) and for the overall data set using Arlequin v3.5.2.2 [43] and DnaSP. The population Phylogenetic and phylogeographical reconstruction genetic structure among AEZs and among haplogroups Different published haplotype sequences of R. appendi- was evaluated by an analysis of molecular variance culatus from South Africa, Kenya, Grande Comore, (AMOVA) performed in Arlequin. The significance of Zambia, Zimbabwe, Uganda and Rwanda were retrieved AMOVA fixation indices was evaluated based on 1,023 from the National Center for Biotechnology Information random permutations. In addition, to assess the level of (NCBI) database (see Additional file 1: Table S1) and genetic distance/differentiation between populations, we were compared with sequences obtained in the present estimated gene flow expressed as of the absolute number study. Phylogenetic reconstruction was performed separ- of migrants (Nm) exchanged among populations, average ately on cox1 and 12S rRNA gene sequences to deter- number of nucleotide differences between populations mine the relationship among populations and possible (K ) and population pairwise genetic differentiation historical dispersal event. To find the evolutionary sub- XY (F ) also in Arlequin [44]. Their significance was tested stitution model that best describe the evolution of cox1 ST using 1000 random permutations. The genetic differenti- and 12S rRNA gene sequences, we performed a hier- ation and population structure statistics were tested archical likelihood ratio test, based on the lowest Bayes- under the hypothesis that different populations are rep- ian information criterion using MEGA v7.0 [49]. The resented by distinct genetic groups or are exchanging nucleotide substitution model selected was then used to Amzati et al. Parasites & Vectors (2018) 11:329 Page 6 of 19 perform a Neighbor-Joining (NJ) and/or Maximum Like- molecular identity (99–100%) with known sequences of lihood (ML) algorithm in MEGA. The stability and R. appendiculatus found in GenBank (Additional file 2: branches support were obtained using 1000 bootstrap Table S2). Haplotype sequences (55 and 23 sequences permutations. Rhipicephalus eversti and Rhipicephalus for cox1and 12S rRNA, respectively) obtained for microplus from this study (GenBank accession numbers each of the six AEZs were deposited and are avail- MF458972 and MF458973 for cox1 and MF479198 and able in GenBank (GenBank accession numbers: MF479199 for 12S RNA genes) and Rhipicephalus tura- MF458895-MF458949 and MF479166-MF479188 for nicus obtained from the GenBank (accession numbers cox1and 12S rRNA genes, respectively). KU880574 and DQ849231 for cox1and 12S rRNA genes, respectively) were used as outgroup taxa. A Phylogenetic relationships and haplotypes distribution Median-Joining (MJ) network was constructed to inves- The overall sequence analysis revealed that cox1had 27 tigate the phylogenetic and ancestral relationship among polymorphic sites, 21 of which were parsimony inform- haplotypes using PopArt Software [50]. ative and 6 were singletons, yielding 22 haplotypes (Additional file 3:Table S3). cox1 amino acid sequences Results did not contain any internal stop codon or indel. Most Morphological and molecular identification of ticks nucleotide mutations were synonymous, except one PCR fragments of the mitochondrial cox1 and 12S rRNA non-synonymous change identified at position 32 of the gene loci were successfully generated from 219 and 126 haplotype CH22 (change from an Alanine to a Threonine). individual ticks, respectively, originating from the three The highest number of cox1 haplotypes was found in AEZs of DRC, the two AEZs of Burundi, one AEZ of DRC AEZ1, while the lowest was observed in Burundi Rwanda and specimens from Tanzania. The 10 se- AEZ3 (Table 2). The 12S rRNA gene provided 10 poly- quences from Tanzania and the sequences previously de- morphic sites, five of which were parsimony informative, scribed in Rwanda were not included in the population defining 9 haplotypes. The 22 cox1 haplotype sequences genetic and structure analyses. Generated nucleotide se- obtained in this study were submitted to GenBank under quences were aligned with the reference haplotype se- accession numbers MF458950-MF458971 and the nine quences retrieved from the GenBank (Additional file 1: 12S rRNA gene haplotypes were deposited under acces- Table S1). The final fragment length obtained for cox1 sion numbers MF479189-MF479197. The phylogenetic re- was 586 bp and 12S rRNA yielded a fragment of 346 bp lationships among cox1haplotypes inferredbya NJ long, with no indels detected in both genes. Their mor- phylogenetic tree (Fig. 2), a ML tree (Fig. 3) and a MJ net- phological identification has been confirmed by the high work (Fig. 4) produced identical topologies and detected Table 2 Rhipicephalus appendiculatus cox1 haplotypes distribution (%) and genetic diversity indices in six agro-ecological zones of the Democratic Republic of Congo, Burundi and Rwanda a,b c Country AEZ n h Haplotype (frequency in %) Haplogroup Genetic diversity indices (%) AB S Hd (SD) K π (SD) (PIS) DRC AEZ1 46 10 CH1(30), CH2(28), CH5(6), CH6(2), CH7(4), CH8(2), CH11(9), CH13(11), 85 15 17 0.82 4.4 (2.2) 0.007 CH16(4), CH17(2) (16) (0.03) (0.002) AEZ2 54 12 CH1(18), CH2(30), CH5(28), CH7(2), CH11(2), CH12(4), CH13(4), 92 8 20 0.81 3.1 (1.6) 0.005 CH18(4), CH19(4), CH20(2), CH21(2), CH22(2) (17) (0.03) (0.001) AEZ3 46 7 CH1(17), CH2(33), CH5(28), CH11(9), CH12(6), CH13(4), CH14(2) 96 4 16 0.79 2.4 (1.3) 0.004 (15) (0.03) (0.001) Burundi AEZ1 26 8 CH1(50), CH2(19), CH3(4), CH4(4), CH5(11), CH6(4), CH7(4), CH10(4) 96 4 18 0.72 2.01 (1.2) 0.003 (3) (0.08) (0.001) AEZ3 17 5 CH1(29), CH2(35), CH5(23), CH8(6), CH9(6) 100 0 4 (2) 0.77 1.1 (0.9) 0.002 (0.06) (0.0003) Rwanda AEZ2 20 8 CH1(30), CH2(20), CH4(5), CH5(5), CH7(20), CH11(5), CH13(10), 70 30 18 0.85 6.3 (3.1) 0.011 CH15(5) (13) (0.05) (0.002) Total 209 22 – 90 10 27 0.81(0.01) 3.4 (1.7) 0.006 (21) (0.0006) Abbreviations: AEZ agro-ecological zones, n number of sequences, h number of haplotypes, S segregation sites, PIS parsimony informative sites, Hd haplotype diversity, SD standard deviation; K mean number of pairwise nucleotide differences, π nucleotide diversity, CH1-22 names of cox1 haplotypes Haplotypes belonging to the haplogroup B are underlined Bold indicates shared haplotypes by all agro-ecological zones A and B are haplogroup names Amzati et al. Parasites & Vectors (2018) 11:329 Page 7 of 19 Fig. 2 Phylogenetic tree of R. appendiculatus cox1 haplotypes. The evolutionary history was inferred by using the neighbor-joining method based on the Tamura 3-parameter model. A discrete Gamma distribution was used to model evolutionary rate differences among sites. Bootstrap values (> 80) are displayed above nodes. CH1-22 are haplotype names. The values in parentheses correspond to the frequency of each haplotype. KU725893 and AF132833 are GenBank accession numbers for R. appendiculatus sequences used as reference haplotypes from Kenya and Zimbabwe, respectively. Two species (R. eversti and R. microplus) obtained in this study and R. turanicus from GenBank (accession number: KU880574) were included as outgroups two distinct clades or lineages strongly-supported by a NJ The distribution of haplotypes presented in Table 2 bootstrap value of 100%. The two lineages diverged at showed that there were shared and private cox1 haplo- least by 12 mutational steps (Additional file 3: Table S3) types confined to restricted AEZs. Three haplotypes but shared a wide range of agro-ecological and geograph- CH1 (27%), CH2 (28%) and CH5 (19%) were shared be- ical conditions in the Great Lakes region. The first lineage, tween all AEZs and were by far the most ubiquitous in named “haplogroup A”, was represented by the most fre- the region, accounting for 74% (154 sequences) of the quent haplotypes (CH1, CH2 and CH5) and comprised in overall dataset (Additional file 3: Table S3). Haplotype total 19 haplotypes consisting of 189 of the 209 sequences CH1 was detected in 13 (50%) of the 26 sequences from analysed (90%), whereas the second lineage “haplogroup lowlands of Burundi (Burundi AEZ1). Two haplotypes B” included three haplotypes (CH7, CH13, CH20) and (CH11 and CH13) defined by 10 and 11 sequences, re- had a total of 20 sequences (10%) (Fig. 2, Additional file 3: spectively, were common in all the AEZs of DRC and in Table S3). Haplogroup B comprised of haplotypes present Rwanda AEZ2. Haplotype CH12 was exclusive to DRC in most AEZs except the highlands region of Burundi AEZ2 and AEZ3. Nine out of the 22 haplotypes were (Burundi AEZ3). These cox1 haplogroups presented a found in single individuals; therefore, they belonged to star-like pattern in the MJ network with less frequent and particular AEZs (Additional file 3: Table S3). Further- single haplotypes connected together to the predominant more, two 12S rRNA haplotypes (12SH1 and 12SH2) or ancestral haplotypes generally by single substitutions were the most abundant, representing 90% of the 126 (Fig. 4), supporting the hypothesis of a recent population sequences analysed. Haplotypes 12SH4 and 12SH5 had expansion. The phylogenetic relationships found in the together 8 (6%) out of the 126 analysed sequences. The cox1 gene were fully congruent with those revealed by the presence of single haplotypes indicates high frequency of ML tree performed on 12S rRNA haplotypes. rare alleles, which suggests a recent population expansion. Amzati et al. Parasites & Vectors (2018) 11:329 Page 8 of 19 Fig. 3 Phylogenetic tree of cox1 haplotypes displaying the relationship between the R. appendiculatus specimens in sub-Saharan African countries. The evolutionary history was inferred by using the maximum likelihood method based on the Tamura 3-parameter model. A discrete Gamma distribution was used to model evolutionary rate differences among sites [5 categories (+G, parameter = 0.39)]. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. Bootstrap scores > 70 are displayed to support nodes. The values in bracket behind haplotype names correspond to the frequency of each haplotype. Haplotype sequences (CH1-22) obtained in the present study are indicated by a black square. Rhipicephalus eversti and R. microplus obtained in this study and R. turanicus (GenBank number: KU880574) were used as outgroups Population genetic diversity from 0.72 in Burundi AEZ1 to 0.85 in Rwanda AEZ2 and Population genetic indices were estimated using cox1nu- the nucleotide diversity (π) values ranged from 0.002 in cleotide sequences and are shown in Table 2. The overall Burundi AEZ3 to 0.011 in Rwanda AEZ2. The average haplotype diversity (Hd) and nucleotide diversity (π)were number of nucleotide pairwise differences (k) was 3.4 for 0.81 and 0.006, respectively. The haplotype diversity ranged the overall dataset, with the highest value observed in Amzati et al. Parasites & Vectors (2018) 11:329 Page 9 of 19 Fig. 4 Median-joining network of the 36 cox1 haplotypes for R. appendiculatus ticks across sub-Saharan African countries. Lines represent mutations and the dot corresponds to a possible intermediate haplotype. Each circle denotes a unique haplotype. Haplotype frequencies are not shown here, but their occurrences across Africa are presented in Table 5 Rwanda AEZ2 (6.3). Altogether, these population genetic (K ), the pairwise genetic distance (F ), and the number XY ST indices showed that the diversity of R. appendiculatus is of migrants (Nm) showed evidence of differentiation quite different among AEZs. Ticks from Rwanda AEZ2 between some R. appendiculatus populations. The popu- were more highly diverse than those from three AEZs of lation pairwise genetic distance (F ) varied from a nega- ST DRC and the two AEZs of Burundi. We also observed an tive and not significant value (F = -0.014, P = 0.87) with ST excess of singleton mutations in Burundi AEZ1 (15 out of infinite value of Nm (between DRC AEZ2 and AEZ3) to 18 polymorphic sites), suggesting a suddenpopulationex- the strongest differentiation (F =0.19, P =0.005)with ST pansion. These data analysed separately by haplogroups low number of migrants (Nm = 2.1) between Burundi (Table 3), showed that the differences of genetic diversity AEZ3 and Rwanda AEZ2. Pairwise F values were found ST among AEZs was largely affected by the frequency distribu- to be low and not significant between DRC AEZ3 and tion of the two cox1 haplogroups (Table 2). Burundi AEZ3 with infinite Nm and between DRC AEZ2 and Burundi AEZ3 associated to a high number of mi- Population structure and ecological differentiation based grants (Nm = 11,875), indicating an excess of gene flow on cox1 haplotypes between these zones. DRC AEZ1 did not differ significantly Analysis of molecular variance (AMOVA) based on cox1 with Burundi AEZ1 (F =0.022, P =0.19) andRwanda ST sequences revealed statistically significant variance among AEZ2 (F = 0.018, P = 0.24). Tick populations from ST the six AEZs analysed together for both combined Rwanda AEZ2 showed a strong genetic differentiation with haplogroups and haplogroup A alone (6%, P <0.001) those from DRC AEZ3 (F =0.18, P <0.001)and DRC ST (Additional file 4: Table S4). The genetic variability among AEZ2 (F =0.14, P = 0.03), indicating reduced or lack of ST individuals within AEZs explained the large majority of gene flow between these populations. These results were molecular variance (94%) for the overall dataset, suggest- confirmed by inter-population nucleotide differences (K ), XY ing an admixture between different populations and the which was significant when populations were significantly coexistence of the two genetically divergent lineages (Fig. 2, differentiated by the pairwise F statistic. Tick populations ST Table 2). This is consistent with moderate genetic struc- were compared by analysing haplogroup A sequences ture of R. appendiculatus in the Great Lakes region. The alone. The findings showed that ticks from Rwanda AEZ2 population differentiation indices were then calculated to and the two AEZs of Burundi were not genetically different. compare the genetic variation observed among AEZs They shared more migrants belonging to haplogroup A (Table 4). Both differentiation statistics based on pairwise (Table 4). In general, the population differentiation was ob- estimates of the Inter-population nucleotide differences served between lowlands and highlands AEZs. Amzati et al. Parasites & Vectors (2018) 11:329 Page 10 of 19 Table 3 Genetic diversity and evolutionary dynamics of the two haplogroups (A and B) identified from cox1 sequences of R. appendiculatus Genetic indices and Haplogroup A Haplogroup Overall data statistics B set DRC Burundi Rwanda Haplogroup AEZ2 A (overall) AEZ1 AEZ2 AEZ3 AEZ1 AEZ3 Diversity indices Number of sequences 39 50 44 25 17 14 189 20 209 Number of polymorphic 79 68 47 18 2 27 sites Number of haplotypes 8 9 6 7 5 6 19 3 22 Haplotype diversity (SD) 0.76 (0.044) 0.78 (0.032) 0.77 (0.033) 0.7 (0.084) 0.77 (0.057) 0.77 (0.089) 0.77 (0.016) 0.56 (0.063) 0.81 (0.01) Nucleotide diversity (SD) 0.002 0.0025 0.002 0.002 0.002 0.002 0.002 0.001 0.006 (0.0003) (0.0003) (0.0003) (0.0004) (0.0003) (0.0006) (0.0001) (0.0002) (0.0006) Neutrality tests Tajima's D (P-value) -0.6 (0.35) -0.79 (0.23) 0.096 (0.6) -1.4 -0.21 (0.45) -1.4 (0.07) -1.5 (0.032)* 0.24 (0.67) -0.93 (0.23) (0.049)* Fu's Fs (P-value) -2.4 (0.08) -2.5 (0.097) -0.1 (0.49) -2.6 -1.1 (0.18) -2.1 -10.4 0.2 (0.49) -3.8 (0.15) (0.027)* (0.041)* (0.001)* Demographic expansion Sum of Squared deviation 0.002 (0.66) 0.004 (0.31) 0.002 (0.49) 0.002 (0.72) 0.02 (0.17) 0.001 (0.9) 0.0008 0.029 0.017 (0.1) (P-value) (0.44) (0.049)* Harpending’s Raggedness 0.06 (0.39) 0.079 (0.21) 0.06 (0.38) 0.05 (0.71) 0.16 (0.12) 0.06 (0.73) 0.057 (0.2) 0.21 (0.1) 0.049 (0.51) index (P-value) Spatial expansion Sum of Squared deviation 0.002 (0.5) 0.004 (0.21) 0.002 (0.4) 0.002 (0.67) 0.02 (0.1) 0.001 (0.9) 0.0008 0.029 0.036 (0.28) (P-value) (0.26) (0.011)* Harpending’s Raggedness 0.064 (0.39) 0.079 (0.21) 0.062 (0.43) 0.05 (0.75) 0.16 (0.12) 0.06 (0.72) 0.057 (0.22) 0.21 (0.12) 0.049 (0.74) index (P-value) Abbreviations: D Tajima’s neutrality statistic, Fs Fu’s neutrality statistic *Values are statistically significant at P < 0.05; significance was determined using 1000 coalescent simulations Demographic and dispersal dynamics of R. appendiculatus haplogroup A, we detected significant evidence of demo- Demographic and spatial dynamics inferred from pairwise graphic and spatial expansion events from the unimodal nucleotide differences revealed bimodal pattern for the mismatch distribution, together with significantly nega- total dataset (Fig. 5a). Tajima’s D and Fu’s Fs were negative tive values for Tajima’s D (D = -1.5, P = 0.032) and Fu’s but not significant (Table 3). However, the sum of square Fs (Fs = -10.4, P = 0.001) statistics. A good fit of sudden deviation (SSD) and raggedness index (RI) for both demo- population expansion was also observed in this hap- graphic and range expansion did not deviate significantly logroup, based on sum of squared deviation values that from that expected under expansion model. The negative were not significant in all the cases: demographic (SSD = values of neutrality tests and the non-significance of SSD 0.0008, P = 0.44) and range (SSD = 0.0008, P = 0.26) ex- and RI indices, suggest a sudden expansion of R. appendi- pansion, with no significant variation of the raggedness culatus populations in the Great Lakes region. The popu- index for both models (Table 3). In contrast, haplogroup lation dynamics history was also analysed separately for B did not display any significant signature of expansion each of the six AEZs (Fig. 6, Additional file 5: Table S5). from the selective neutrality tests. In addition, the ob- Most AEZs showed a bimodal mismatch distribution, served mismatch pattern for this haplogroup deviated except in Burundi AEZ3 where the mismatch pattern was significantly from that expected under population expan- unimodal (Fig. 6e). The observed bimodal pattern sug- sion scenario (SSD = 0.029, P = 0.049 and P = 0.011 for gested population subdivision as shown by the existence demographic and spatial expansion, respectively), imply- of two well resolved haplogroups (Fig. 2). ing that haplogroup B did not experience any historical The population dynamics were then inferred for each population expansion. This group is characterised either haplogroup separately (Table 3). A unimodal distribution by a constant population size (demographic equilibrium) was observed in both haplogroups (Fig. 5b, c). For the or had experienced a weak signal of population Amzati et al. Parasites & Vectors (2018) 11:329 Page 11 of 19 Table 4 Population genetic statistics for pairwise comparison of different populations of R. appendiculatus from sequences of cox1 gene. Values in parentheses represent the P-value statistics Population 1 Population 2 Haplogroup A and B Haplogroup A K F Nm K F Nm XY ST XY ST DRC AEZ1 DRC AEZ2 3.8 (0.044)* 0.036 (0.044)* 13.2 1.5 (< 0.001)* 0.11 (< 0.001)* 4 DRC AEZ1 DRC AEZ3 3.6 (0.032)* 0.057 (0.032)* 8.3 1.5 (0.001)* 0.1 (0.001)* 4.3 DRC AEZ1 Burundi AEZ1 3.3 (0.15) 0.022 (0.19) 22.2 1.2 (0.31) 0.002 (0.33) 243.6 DRC AEZ1 Burundi AEZ3 3 (0.059) 0.056 (0.087) 8.4 1.3 (0.087) 0.043 (0.082) 11.2 DRC AEZ1 Rwanda AEZ2 5.4 (0.3) 0.018 (0.24) 26.8 1.3 (0.78) -0.026 (0.79) ∞ DRC AEZ2 DRC AEZ3 2.7 (0.85) -0.014 (0.87) ∞ 1.4 (0.93) -0.016 (0.93) ∞ DRC AEZ2 Burundi AEZ1 2.7 (0.003)* 0.060 (0.011)* 7.9 1.5(< 0.001)* 0.13 (< 0.001)* 3.2 DRC AEZ2 Burundi AEZ3 2.1 (0.23) 0.00004 (0.43) 11876 1.3 (0.46) -0.012 (0.54) ∞ DRC AEZ2 Rwanda AEZ2 5.2 (0.003)* 0.14 (0.003)* 3.1 1.6 (0.006)* 0.11 (0.004)* 3.8 DRC AEZ3 Burundi AEZ1 2.4 (< 0.001)* 0.076 (0.006)* 6.1 1.5(< 0.001)* 0.14 (< 0.001)* 3.2 DRC AEZ3 Burundi AEZ3 1.8 (0.32) -0.005 (0.43) ∞ 1.3 (0.39) -0.009 (0.46) ∞ DRC AEZ3 Rwanda AEZ2 5.1 (< 0.001)* 0.18 (< 0.001)* 2.2 1.6 (0.009)* 0.11 (0.011)* 3.9 Burundi AEZ1 Burundi AEZ3 1.6 (0.074) 0.034 (0.08) 12.2 1.2 (0.077) 0.063 (0.073) 7.4 Burundi AEZ1 Rwanda AEZ2 4.8 (0.030)* 0.14 (0.03)* 3.1 1.2 (0.97) -0.039 (0.98) ∞ Burundi AEZ3 Rwanda AEZ2 4.6 (0.005)* 0.19 (0.005)* 2.1 1.3 (0.16) 0.044 (0.14) 10.9 Abbreviations: K average number of nucleotide differences between populations, F pairwise genetic distance F-statistic based on nucleotide sequences (Right’s XY ST fixation index), Nm number of migrants between populations *Values are statistically significant at P < 0.05; significance was determined using 1000 coalescent simulations. bottleneck that reduced its diversity. When analysing the each of the six AEZs and for Tanzania (Table 5). demographic dynamics for samples belonging to haplogroup Twenty-eight haplotypes had been already described in dif- A in each AEZ (Table 3), population expansion signal was ferent African countries, and eight haplotypes (CH3, CH9, confirmed in all AEZs by mismatch analyses exhibiting uni- CH10, CH14, CH16, CH18, CH20 and CH33) were newly modal distribution (Additional file 6:FigureS1).For thesix described in the present study. Most of these new haplo- AEZs, none of the statistical comparisons between the ob- types were less abundant (Additional file 3:Table S3)and served and the expected distributions rejected the sudden diverged from the common ancestral haplotypes generally and range expansion assumptions based on the raggedness by only one substitution (Fig. 4). Haplotypes CH1, CH7, index and the sum of squared deviation. The neutrality indi- CH11 were the most ubiquitous and shared wide geo- ces were generally negative but not significant, except in graphical distribution in affected African countries. CH1 Burundi AEZ1 where both Tajima’s D (D =-1.4, P = 0.049) haplotype was shared by Kenya, Eastern Zambia, DRC and Fu’s Fs (Fs =-2.6, P = 0.027) statistics showed significant (AEZ1, 2 and 3), Burundi (AEZ1 and 3), Tanzania and negative values and in Rwanda AEZ2 where Fu’s Fs was Rwanda (AEZ2 and GenBank sequences), while haplotype negative (Fs = -2.1) and significant (P = 0.041). According to CH7 was reported in Kenya, South Africa, Zimbabwe, the mismatch distribution and negative values for neutrality Grande Comore, Eastern and southern provinces of tests, the hypothesis of population growth and spatial expan- Zambia, DRC (AEZ1, 2 and 3), Burundi (AEZ1) and sion models could not be rejected in the six AEZs, which Rwanda (AEZ2). Haplotype CH11 was present in Kenya, was consistent with a model of sudden expansion for each Rwanda (AEZ2 and GenBank sequences), Comoros and population subdivision. DRC (AEZ1, 2 and 3). Eighteen haplotypes mostly with unique sequences were found to be restricted to Kenya Phylogeographical structure and have not been reported in any other country. In the To study the phylogeographical structure of R. appendicu- same way, haplotype CH23 was particular to Uganda. This latus in Africa, the haplotype sequences found in this study country did not share any haplotype with other countries. along with those retrieved from GenBank (Additional file 1: The ML phylogenetic tree reconstructed using the 105 se- Table S1) were used to reconstruct the phylogenetic using quences indicated that globally, the African tick R. appen- ML tree and the MJ network methods based on cox1and diculatus is consistently clustered into two groups 12S rRNA genes. Thirty-six cox1 haplotypes were identified (haplogroups A and B) well-supported by a ML bootstrap in 105 sequences including 50 haplotypes from GenBank value of 100% (Fig. 3). Our 19 haplotypes that had been and the 55 haplotype sequences obtained in our study for described for haplogroup A (Table 2,Additional file 3: Amzati et al. Parasites & Vectors (2018) 11:329 Page 12 of 19 Fig. 5 cox1 mismatch distribution pattern for two haplogroups of R. appendiculatus. a shows the overall mismatch distribution pattern for all AEZs, b and c show the mismatch distribution of nucleotide sequences in haplogroups A and B, respectively. The x-axis shows the number of pairwise differences between pairs of haplotype sequences and the y-axis shows their frequencies (in %). The observed frequencies are represented by solid histograms and the simulated mismatch distributions expected under demographic expansion (solid black line) and under spatial expansion (dotted black line). Simulated curves under range and demographic expansion have same pattern in these figures, they overlapped Table S3) formed one clade with 19 haplotypes from Kenya (19/29), all the 3 haplotypes from Rwanda (3/3), six haplotypes from Eastern province of Zambia (6/7), all the five haplotypes from Tanzania, whereas our three haplo- types of haplogroup B were clustered with 10 haplotypes from Kenya, the single haplotype from Grande Comore, all haplotypes from South Africa (3/3), Zimbabwe (3/3), Uganda (2/2), Southern province of Zambia (2/2), and one haplotype from Eastern province of Zambia. In addition to ML tree, the MJ network also revealed that R. appendi- culatus is divided into two main groups in Africa, sepa- rated by seven mutational steps (Fig. 4). Similar findings were confirmed by 12S rRNA gene. Our 23 12S rRNA individual haplotypes from each of the six AEZs were analysed together with the 29 sequences ob- tained from GenBank (Additional file 1: Table S1). Fourteen haplotypes were observed, two most common (12SH1 and 12SH5) and 12 minors (defined mostly by one sequence or restricted to particular country) (Additional file 7:Table S6). Haplotype 12SH1 was common in DRC, Burundi, Rwanda, Kenya and Eastern province of Zambia, while haplotype 12SH5 was present in DRC, Rwanda, Zimbabwe, Comoros, South Africa, Eastern Zambia and Kenya. Six new haplo- types were not found outside the Great lakes region (12SH3, 12SH4 and 12SH6-H9). The NJ phylogenetic resolved these 12S rRNA haplotypes into two haplogroups (haplogroup A and B) supported by 97% bootstrap value (Additional file 8: Figure S2). Their pattern was identical to that observed from cox1 haplogroups. However, the Ugandan haplotype (12SH10: GenBank AF150028) was clustered in haplogroup A, showing that the haplogroup A is also present in Uganda. Unfortunately, we did not find its corresponding cox1 sequence. Discussion This study analysed the intraspecific variation of mitochon- drial DNA to better understand how factors such as agro-ecological zones and anthropogenic movements of cat- tle may affect population genetic structure and population expansion history of the tick R. appendiculatus,the main vector of T. parva in sub-Saharan Africa. We expected evi- dence of R. appendiculatus population expansion, gene flow Amzati et al. Parasites & Vectors (2018) 11:329 Page 13 of 19 Fig. 6 cox1 mismatch distribution pattern for six populations of R. appendiculatus. a, b and c show the mismatch distribution pattern for R. appendiculatus from DRC AEZs (AEZ1, 2 and 3, respectively); d and e represent the mismatch pattern of ticks from Burundi AEZ1 and AEZ3, respectively; f depicts the mismatch distribution of ticks from Rwanda AEZ2. The x-axis shows the number of pairwise differences between pairs of haplotype sequences and the y-axis shows their frequencies in %. The observed frequencies are represented by solid histograms. Black full line represents the expected distribution under sudden expansion model, and dotted line represents the distribution simulated under spatial expansion model. Simulated curves under spatial and demographic expansion have same pattern in (d), and they overlapped and different colonization patterns of tick lineages, facilitated groups [29], corresponding to our haplogroups A and B, by the reported cattle mobility in the Great Lakes region. respectively. This genetic grouping fitted well with the phenotypical, physiological and ecological variation stud- Two R. appendiculatus lineages that are more variable in ies which have previously distinguished two major sub- lowlands than highlands occur in the Great Lakes region populations of R. appendiculatus in Africa [32–37, 51]. The 22 haplotypes identified by DNA polymorphic ana- These phenotypic and physiological variations are largely lysis of cox1 gene locus were clustered into two associated to agro-ecological and geographical subdivi- well-defined major groups, named haplogroup A (the sions. Tropical areas with prolonged and marked dry most frequent) and haplogroup B. Similar grouping were seasons are more suitable for larger sized ticks express- obtained with 12S rRNA analyses. The two haplogroups ing high intensity of diapause and displaying univoltine identified in the present study have been previously de- phenology, corresponding to “southern African group” scribed as “east African” and “southern African” genetic or haplogroup B. Equatorial areas with bimodal or Amzati et al. Parasites & Vectors (2018) 11:329 Page 14 of 19 Table 5 Rhipicephalus appendiculatus cox1 haplotypes and their distribution among agro-ecological zones of the Great lakes region and other sub-Saharan African countries Haplotype Haplotypes from GenBank: Country (original haplotype name and Present study Haplogroup GenBank number) a e CH1 Kenya (H2: KU725891, H4: KU725893, H6: KU725895) , Rwanda Burundi (AEZ1, AEZ3), DRC (AEZ1, AEZ2, AEZ3), A f e g (H3:DQ901360) Zambia-east (H4: KX276942 , H5: DQ859266 , Tanzania (TZ18, TZ10, TZ08, TZ20), Rwanda f g H3:DQ901361 , H2: DQ859265 ) (AEZ2) CH3 – Burundi AEZ1 A CH4 Rwanda (H6: DQ901362) Burundi AEZ1, Rwanda (AEZ2) A b e CH6 Kenya (H5: KU725894) Burundi (AEZ1, AEZ3), DRC (AEZ1, AEZ2) A c e e e CH7 Kenya (H1: KU725890) , South Africa (H1: KX276939 , H1: KX276940 , H1: Burundi AEZ1, DRC (AEZ1, AEZ2, AEZ3), Rwanda B f g e DQ901356 ), Zambia-east (H1: DQ859261) , Zambia-south (H1:KX276943 , (AEZ2) g h i e H1: DQ859262) , Zimbabwe (AF132833 , KC503257 , H1: KX276944 ), Grande Comore (H1) CH9 – Burundi AEZ3 A CH10 – Burundi AEZ1 A d e f j CH11 Kenya (H3: KU725892) , Rwanda (H3: DQ901363) , Grande Comore (H3) DRC (AEZ1, AEZ2, AEZ3), Rwanda (AEZ2) A CH12 Kenya (H11: KU725900) DRC (AEZ2, AEZ3) A CH14 – DRC AEZ3 A CH16 – DRC AEZ1 A CH18 – DRC AEZ2 A CH20 – DRC AEZ2 B CH23 Uganda (H8: KX276941, KU725897) – B CH24 Kenya (H14: KU725903) – B CH25 Kenya (H27: KU725916) – B f,j e CH26 Grande Comore (H2: DQ901357) , Kenya (H7: KU725896, H13: KU725902) – B CH27 Kenya (H16: KU725905) – B e f f CH28 Kenya (H21: KU725910 , H9: DQ901359 , H9: DQ901358 ) – B CH29 Kenya (H28: KU725917) – B CH30 Kenya (H9: KU725898) – A CH31 Kenya (H17: KU725906 – A CH32 Kenya (H24: KU725913) – A CH33 – Tanzania (TZ13) A CH34 Kenya (H22: KU725911) – A CH35 Kenya (H23: KU725912) – A CH36 Kenya (H10: KU725899) – A CH37 Kenya (H15:KU725904) – A CH38 Kenya (H20: KU725909) – A CH39 Kenya (H19: KU725908) – A CH40 Kenya (H26: KU725915) – A CH41 Kenya (H25: KU725914 – A CH42 Kenya (H18: KU725907) – A CH43 Kenya (H12: KU725901) – A CH44 Zambia-east (H3: DQ859263) – A CH45 Zambia-east (H4: DQ859264) – A CH1: variants CH1, CH2, CH5, CH15, CH17 and CH21 (Additional file 3: Table S3) CH6: variants CH6, CH8 and CH22 (Additional file 3: Table S3) CH7: variants CH7 and CH13; d CH11: CH11 and CH19 (Additional file 3: Table S3) Haplotypes underlined are exclusive to the Great Lakes region. Similar data for 12S rRNA are detailed in Additional file 7: Table S6 e f g h i j [31]; [29]; [52]; [68]; [69]; [14] Amzati et al. Parasites & Vectors (2018) 11:329 Page 15 of 19 continuous rainfall rather harbour smaller ticks without tick R. appendiculatus that took place in the Great Lakes diapause with bivoltine or continuous phenology, corre- region. A strong evidence of recent spatial and demo- sponding to “east African group” or haplogroup A [37, 52]. graphic expansion was observed for lineage A in all AEZs The highest genetic variability was observed in low- included in the study. Analyses of cox1 sequences revealed lands, whereas a relatively lower diversity was observed relatively high haplotype diversity contrasted with low nu- in midlands and highlands. The high genetic diversity in cleotide diversity values for each population, suggesting a lowlands can be explained by the strong presence of the sudden population expansion [54, 55]. This result, together two lineages A and B observed in these AEZs. The coex- with negative values of Tajima’s D and Fu’s Fs, the star-like istence of these lineages could originate from the disper- radiation, the unimodal mismatch pattern and sal of the tick through livestock movement between non-significant RI and SSD statistics, further support the re- AEZs [53], associated to the suitability of semi-arid cli- cent evolutionary history andsuddenpopulationgrowth mate for lineage B expressing obligatory diapause [37]. experienced by lineage A [56–58]. We did not estimate the expansion time, but we hypothesize that the expansion was Moderate genetic structure of R. appendiculatus between recent and not sufficient to increase the nucleotide diver- lowlands and highlands sity, probably because of recent coalescence, while rapid Population genetic analyses of cox1 gene variation in R. population expansion following a selective sweep (bottle- appendiculatus revealed low to moderate genetic differen- neck or genetic drift) could have accumulated new muta- tiation values and high gene flow rates among AEZs. The tions that sufficiently increased the haplotype diversity [45, two identified R. appendiculatus lineages were sympatric 59, 60]. This could explain the excess of singletons poly- in the Great Lakes region, although lineage A was the morphic sites and rare haplotypes that diverge from most abundant and widely distributed in all AEZs and ancestral haplotypes by only 1–2mutationalsteps [45, 61]. lineage B was particular confined to lowlands, where the However, the strong bimodal mismatch pattern observed in climate is tropical dry, more arid with lower annual rain- Rwanda AEZ2 and DRC AEZ1 suggests the coexistence of fall and longer dry season than in highlands (short dry sea- R. appendiculatus lineages after recent colonization events son with abundant annual rainfall). These climate or exchanging migrants [45, 56, 62]. The scenario of conditions in lowlands are quite similar to the described lineage B contrasts with that of lineage A. Analyses indi- ecological zones of lineage B in southern Africa [29, 52]. cate, for lineage B, an equilibrium situation that is charac- In addition, altitudinal gradient seems to be a key factor terised by a weak signal of recent bottleneck and no that shapes the distribution pattern and the establishment evidence of population expansion. It was also less diverse ability of lineage B. Its frequency decreases with increasing than haplogroup A, indicating that haplogroup A is ex- altitude. High degree of genetic similarity was observed periencing population expansion independently of hap- between the lowlands of DRC and the low plateau of logroup B and it has been established longer in the Great Rwanda and between the highlands of DRC and Burundi, Lakes region, while haplogroup B was probably which are geographically distant from each other. The introduced more recently and established a founder popu- most likely explanation for this is that the spatial pattern lation. There are three possible explanations of the equi- of R. appendiculatus lineages is not only driven by geo- librium observed in haplogroup B: (i) only few haplotypes graphical separation as described in previous studies [52], were recently introduced; (ii) only the few identified hap- but also related to their ecological preferences, as ob- lotypes persisted after a bottleneck; or (iii) haplogroup B is served by the significant genetic differentiation among experiencing an initial selection sweep which has reduced lowlands and highlands AEZs. On the other hand, adja- the number of rare haplotypes and singleton mutations cent AEZs shared more migrants, especially of lineage A, [56]. When we analysed R. appendiculatus cox1sequences facilitated by short-distance seasonal movement of cattle available in Africa, the star radiation of the MJ network in [9], which may have reduced the geographical structuring each group suggests that the two lineages went through a of the tick [11, 18]. Analysis of molecular variance demographic expansion and evolve independently of each (AMOVA) confirmed these findings showing that the vari- other with limited gene exchange. Unfortunately, we were ance explained by divergence between the six AEZs was not able to test the hypothesis of crossbreeding events be- lower (6%), while the largest fraction of genetic variation tween the two lineages described here, because of the ma- was observed among individuals within AEZs (94%). ternal inheritance of mitochondrial genes used in the present study. More studies, such as extensive biological Rhipicephalus appendiculatus lineage A has undergone characterisation, crossbreeding experiments and the use sudden demographic and range expansion in the Great of biparental inheritance markers are necessary to investi- Lakes region gate the panmixia of the two lineages and the genetic The demographic and spatial dynamics were analysed using mechanisms driving their establishment and correspond- multiple algorithms, to elucidate colonization events of the ing phenotypic variations in changed environments. Amzati et al. Parasites & Vectors (2018) 11:329 Page 16 of 19 The fact that R. appendiculatus has undergone spatial behaviour for southern African ticks, which allow them expansion was in accordance with the expectation that to survive hot and dry ecological conditions [33, 37]. We long and short distance movement of cattle are key factors hypothesise that these characteristics make the develop- of spreading ticks. The processes responsible for that evo- ment of lineage B slower than lineage A in sympatric lutionary pattern may have resulted, not only from range areas, giving an evolutionary advantage to the latter. expanding of the previously established haplotypes to This could also reduce the abundance of lineage B delay- proximate AEZs, but also from the recolonization events ing its oviposition during unfavourable conditions [37]. of ticks from other regions and countries [53, 63]. In The processes that took place to divide the two groups addition to cattle movement, the population expansion need to be further investigated. However, we propose and establishment of novel haplotypes in previously un- that they could have diverged following genetic drift due occupied areas could have been driven by recent environ- to founder events of natural geographical barrier that mental and climatic changes, affecting vector-borne may result in reproductive isolation. It is demonstrated diseases landscape over recent decades [64–66]. that biogeographical break again host migration reduces gene exchange and could dictate the spatial and repro- Sympatric and allopatric ecological zones of R. ductive separation within a species [60]. appendiculatus lineages in Africa Phylogenetic trees produced two main genetic subpopu- Conclusions lations of R. appendiculatus that have a wide distribution This study provided new insights into the genetic struc- range in Africa, with large divergence in behavioural dia- ture and colonization events of R. appendiculatus in the pause [36, 37], spatial variation in body size [35, 51], Great Lakes region and over its distribution range in ecology and phenology [32–34]. Initially, the two line- sub-Saharan Africa. Our results highlighted the occur- ages were considered as “east African” and “southern Af- rence of two major genetic lineages (A and B) in the Great rican” genetic groups [29]. To date, they are sympatric Lakes region. The two lineages are not spatially structured in most eastern and central African ecological zones. For in the study region and they differ in their colonization instance, previous studies did not reveal the presence of histories and pattern. Lineage B, not previously reported, haplogroup B in Rwanda [29]. This could be an indica- was probably introduced recently in the region and its oc- tion of recent introduction of the tick in the Great Lakes currence decreases with increased altitude, whereas region. The MJ network further elucidate that the initial lineage A, widely distributed, has been longer established population of haplogroup B could have come from an and subjected to sudden demographic and spatial expan- ancestral female of haplotype CH7, which is the most sion most likely related to short and long-distance cattle prevalent haplotype of this group in areas where it oc- movement. Rhipicephalus appendiculatus ticks are more curs. Consequently, two different eco-genetic zones are diverse in lowlands than highlands with moderate genetic shaped in Africa, the sympatric zone where the two line- differentiation between the two ecosystems, while more ages are found, which covers central and eastern Africa, genetic similitude is found in zones with same and allopatric zone in southern Africa where the two lin- agro-ecological profiles, in spite of their geographical dis- eages have clear geographical and ecological separation. tance. The genetic distribution of R. appendiculatus sug- For instance, in Zambia, lineage B is present in the south gests two different eco-genetic zones in Africa, the (long dry season) and lineage A in the east of the coun- sympatric zone (central and eastern Africa) where the two try (shorter dry season) [29]. In Grande Comore, lineage lineages coexist and the allopatric zones (southern Africa) B has established a stable population, while lineage A where they have clear geographical divergence. The range was found on imported cattle [14]. In areas where the expansion pattern of lineages and the genetic admixture two lineages are sympatric, their respective abundances of R. appendiculatus populations observed in the Great differ, mainly driven by their divergence in ecological Lakes region can strongly affect the epidemiological dy- preferences and plasticity [33]. The sympatric relation- namics of ECF. This could partially explain the endemic ship is in agreement with the observation made by Berk- instability and occasional epidemics due to the introduc- vens et al. [67], where an east African stock from Kenya tion and temporal subsistence of infected ticks mostly in expressed diapause, contrasting with the result obtained fringes areas of lowlands. by Madder et al. [36], where another stock from the same region was unable to enter diapause. These evi- Additional files dences show that lineage B has a greater invasive ability into new habitats and better fit wide range of tropical Additional file 1: Table S1. Rhipicephalus appendiculatus cox1 and 12S rRNA haplotype sequences retrieved from GenBank. (DOCX 16 kb) and equatorial conditions, while lineage A is particularly Additional file 2: Table S2. cox1 and 12S rRNA BLAST results for species confined to equatorial conditions. This could be ex- identification and confirmation. (DOCX 16 kb) plained by the larger body size and obligatory diapause Amzati et al. Parasites & Vectors (2018) 11:329 Page 17 of 19 GenBank under accession numbers MF458950-MF458971 and the GenBank Additional file 3: Table S3. Polymorphism in the 22 haplotypes of the numbers for the nine 12S rRNA haplotypes are MF479189-MF479197. cox1 gene fragment of R. appendiculatus. (DOCX 21 kb) Additional file 4: Table S4. Population genetic structure inferred by Authors’ contributions analysis of molecular variance (AMOVA) based on cox1 sequences of R. GSA collected and identified ticks, performed experiments, analysed data appendiculatus from different agro-ecological zones. (DOCX 13 kb) and wrote the manuscript. TM, NK and GSA conceived the study. TM, NK, RP, JBM, AD and MM participated in the study design, supervised the research Additional file 5: Table S5. Evolutionary neutrality, demographic and and critically revised the manuscript. RP and EGK participated in laboratory spatial history of mitochondrial cox1 gene. (DOCX 15 kb) experiments and advised in data analyses. JBM supervised tick collections. All Additional file 6: Figure S1. cox1 mismatch distribution pattern for R. authors read and approved the final manuscript. appendiculatus haplogroup A in different agro-ecological zones. (DOCX 193 kb) Additional file 7: Table S6. Rhipicephalus appendiculatus 12S rRNA Ethics approval and consent to participate haplotypes and their distribution among agro-ecological zones of the The research was carried out in accordance with ethical guidelines and Great lakes region and other sub-Saharan African countries. (DOCX 15 kb) regulations of the Université Evangélique en Afrique (UEA/Bukavu-DR Congo). The specific field protocols were approved by the UEA (Reference Additional file 8: Figure S2. Neighbor-joining tree of 12S haplotype number: UEA/SGAC/KM/020/2015) and permitted by the Provincial sequences for R. appendiculatus across African countries. (DOCX 18 kb) Inspection of Agriculture, Livestock and Fisheries (Reference number: 55.00/ 0026/IPAPEL/SK/2015). Farm managers and owners were informed about the Abbreviations survey and gave approval for tick collection on their cattle. This research did AMOVA: Analysis of molecular variance; BecA: Biosciences eastern and central not involve species at risk of extinction. Africa; BLAST: Basic local alignment search tool; cox1: Cytochrome c oxidase subunit 1; D: Tajima’s neutrality test; ECF: East Coast fever; Fs:Fu’s neutrality Competing interests tests; F : Right’s fixation index; h: Number of haplotypes; Hd: Haplotype ST The authors declare that they have no competing interests. diversity; ILRI: International Livestock Research Institute; ITS2: Internal transcribed spacer 2; K: Mean number of pairwise nucleotide differences Publisher’sNote within population; K : Average number of nucleotide differences between XY Springer Nature remains neutral with regard to jurisdictional claims in published populations; ML: Maximum likelihood; NCBI: National Center for maps and institutional affiliations. Biotechnology Information; NJ: Neighbor-joining; Nm: Number of migrants; PCR: Polymerase chain reaction; PIS: Parsimony informative sites; Author details RI: Harpending’s raggedness index; rRNA: ribosomal RNA; S: Segregating sites; Unit of Integrated Veterinary Research, Department of Veterinary Medicine, SD: Standard deviation; SSD: Sum of squares deviation; π: nucleotide diversity Faculty of Sciences, Namur Research Institute for Life Sciences (NARILIS), University of Namur (UNamur), Rue de Bruxelles 61, 5000 Namur, Belgium. Acknowledgments Research Unit of Veterinary Epidemiology and Biostatistics, Faculty of We would like to thank Professor Paul Gwakisa from the Sokoine University Agricultural and Environmental Sciences, Université Evangélique en Afrique, of Agriculture (SUA) in Tanzania for providing tick specimens from Tanzania. P.O. Box 3323, Bukavu, Democratic Republic of the Congo. Biosciences We thank Alphone Bisusa from the Centre de Recherche en Sciences eastern and central Africa - International Livestock Research Institute Naturelles (CRSN/Lwiro) and Amzati Saidi Ngton (Université Evangélique en (BecA-ILRI) hub, P.O. Box 30709-00100, Nairobi, Kenya. Department of Afrique) for technical assistance with sample collection and tick Biochemistry, School of Medicine, University of Nairobi, P.O. Box identification. We are also grateful to the technical staff of BecA-ILRI for facili- 30197-00100, Nairobi, Kenya. Present address: Centre for Tropical Livestock tating laboratory work, especially Stephen Mwaura (tick unit-ILRI) for assist- Genetics and Health (CTLGH), The University of Edinburgh, Easter Bush, ance with the confirmation of morphological identification of closely related Midlothian, Scotland EH25 9RG, UK. Department of Veterinary Tropical tick species. We thank Dr Marie Cariou (Research Unit in Environmental and Diseases, Faculty of Veterinary Science, University of Pretoria, P/Bag X04, Evolutionary Biology, University of Namur) for advising on data analysis and Onderstepoort 0110, South Africa. interpretation. We are also grateful to Professor Gustave Mushagalusa (Uni- versité Evangélique en Afrique) for mobilizing partnerships and collabora- Received: 11 January 2018 Accepted: 16 May 2018 tions to the “Theileria” project in DRC. Funding References The laboratory aspects of this project were fully supported by the BecA-ILRI 1. Nene V, Kiara H, Lacasta A, Pelle R, Svitek N, Steinaa L. The biology of Hub through the Africa Biosciences Challenge Fund (ABCF) program. The Theileria parva and control of East Coast fever - Current status and future ABCF Program is funded by the Australian Department for Foreign Affairs trends. Ticks Tick Borne Dis. 2016;7:549–64. and Trade (DFAT) through the BecA-CSIRO partnership; the Syngenta Foun- 2. Kalume MK, Saegerman C, Mbahikyavolo DK, Makumyaviri AM, Marcotty T, dation for Sustainable Agriculture (SFSA); the Bill & Melinda Gates Foundation Madder M, et al. Identification of hard ticks (Acari: Ixodidae) and (BMGF); the UK Department for International Development (DFID) and; the seroprevalence to Theileria parva in cattle raised in North Kivu Province, Swedish International Development Cooperation Agency (Sida). Sample col- Democratic Republic of Congo. Parasitol Res. 2013;112:789–97. lection, field equipment and scientific workshop organised in DRC were sup- 3. Bazarusanga T, Geysen D, Vercruysse J, Madder M. An update on the ported through the “Theileria” project co-funded to the Université ecological distribution of ixodid ticks infesting cattle in Rwanda: Evangélique en Afrique (UEA) by the Agence Universitaire de la Francopho- countrywide cross-sectional survey in the wet and the dry season. Exp Appl nie (AUF) and the Communauté Economique des Pays des Grands Lacs Acarol. 2007;43:279–91. (CEPGL). The International foundation of Science (IFS) supported in collabor- 4. Kaiser M, Sutherst R, Bourne A, Gorissen L, Floyd R. Population dynamics of ation with BecA-ILRI hub the individual scholarship awarded to GSA for field ticks on Ankole cattle in five ecological zones in Burundi and strategies for work and part of field equipment to the “Theileria” project. The study was their control. Prev Vet Med. 1988;6:199–222. also supported by the University of Namur through the UNamur-CERUNA in- 5. Olwoch J, Reyers B, van Jaarsveld A. Host-parasite distribution patterns stitutional PhD grant awarded to GSA for manuscript write-up in Belgium. under simulated climate: implications for tick-borne diseases. Int J Climatol. 2009;29:993–1000. Availability of data and materials 6. Perry B, Lessard P, Norval R, Kundert K, Kruska R. Climate, vegetation The data supporting the findings of this study are included within the article and the distribution of Rhipicephalus appendiculatus in Africa. Parasitol and its additional files. Representative sequences of R. appendiculatus from Today. 1990;6:100–4. each agro-ecological zone (AEZ) are available in the GenBank database under 7. Mararo SB. Pouvoirs, élevage bovin et la question foncière au Nord-Kivu. In: accession numbers: MF458895-MF458949 and MF479166-MF479188 for cox1 Marysse S, Reyntjens F, editors. L'Afrique des Grands Lacs: annuaire 2000– and 12S rRNA genes, respectively. The 22 cox1 haplotypes were submitted to 2001. Paris: L'Harmattan; 2001. p. 219–50. Amzati et al. Parasites & Vectors (2018) 11:329 Page 18 of 19 8. De Failly D. L’économie du Sud-Kivu 1990–2000: mutations profondes 33. Speybroeck N, Madder M, Van Den Bossche P, Mtambo J, Berkvens N, Chaka cachées par une panne. In: Marysse S, Reyntjens F, editors. L’Afrique des G, et al. Distribution and phenology of ixodid ticks in southern Zambia. Med Grands Lacs: annuaire 1999–2000. Paris: L'Harmattan; 2000. p. 163–92. Vet Entomol. 2002;16:430–41. 9. Verweijen J, Brabant J. Cows and guns. Cattle-related conflict and 34. Leta S, De Clercq EM, Madder M. High-resolution predictive mapping for armed violence in Fizi and Itombwe, eastern DR Congo. J Mod Afr Rhipicephalus appendiculatus (Acari: Ixodidae) in the Horn of Africa. Exp Appl Stud. 2017;55:1–27. Acarol. 2013;60:531–42. 10. Vlassenroot K, Huggins C. Land, migration and conflict ineasternDRC. In: 35. Chaka G, Billiouw M, Geysen DM, Berkvens DL. Spatial and temporal Huggins C, Clover J, editors. From the ground up: land rights, conflict and peace variation in Rhipicephalus appendiculatus size in eastern Zambia. Trop Med in sub-Saharan Africa. Pretoria: Institute for Security Studies; 2005. p. 115–94. Int Health. 1999;4:A43–8. 11. Boulinier T, Kada S, Ponchon A, Dupraz M, Dietrich M, Gamble A, et al. 36. Madder M, Speybroeck N, Brandt J, Berkvens D. Diapause induction in Migration, prospecting, dispersal? What host movement matters for adults of three Rhipicephalus appendiculatus stocks. Exp Appl Acarol. 1999; infectious agent circulation? Integr Comp Biol. 2016;56:330–42. 23:961–8. 12. Maze-Guilmo E, Blanchet S, McCoy KD, Loot G. Host dispersal as the driver 37. Madder M, Speybroeck N, Brandt J, Tirry L, Hodek I, Berkvens D. Geographic of parasite genetic structure: a paradigm lost? Ecol Lett. 2016;19:336–47. variation in diapause response of adult Rhipicephalus appendiculatus ticks. 13. De Deken R, Martin V, Saido A, Madder M, Brandt J, Geysen D. An outbreak Exp Appl Acarol. 2002;27:209–21. of East Coast fever on the Comoros: a consequence of the import of 38. Bazarusanga T, Vercruysse J, Marcotty T, Geysen D. Epidemiological studies immunised cattle from Tanzania? Vet Parasitol. 2007;143:245–53. on theileriosis and the dynamics of Theileria parva infections in Rwanda. Vet 14. Yssouf A, Lagadec E, Bakari A, Foray C, Stachurski F, Cardinale E, et al. Parasitol. 2007;143:214–21. Colonization of Grande Comore Island by a lineage of Rhipicephalus 39. Walker AR, Bouattour A, Camicas J-L, Estrada-Pena A, Horak IG, Latif AA, appendiculatus ticks. Parasit Vectors. 2011;4:38. et al. Ticks of domestic animals in Africa: a guide to identification of species. Edinburgh, Scotland: Bioscience Reports; 2003. 15. Fevre EM, Bronsvoort BM, Hamilton KA, Cleaveland S. Animal movements and the spread of infectious diseases. Trends Microbiol. 2006;14:125–31. 40. Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R. DNA primers for 16. Barre N, Uilenberg G. Spread of parasites transported with their hosts: case amplification of mitochondrial cytochrome c oxidase subunit I from diverse study of two species of cattle tick. Rev Sci Tech. 2010;29:135–47. metazoan invertebrates. Mol Mar Biol Biotechnol. 1994;3:294–9. 17. Estrada-Peña A, Salman M. Current limitations in the control and spread of 41. Simon C, Frati F, Beckenbach A, Crespi B, Liu H, Flook P. Evolution, ticks that affect livestock: a review. Agriculture. 2013;3:221–35. weighting, and phylogenetic utility of mitochondrial gene sequences and a 18. Leo SS, Gonzalez A, Millien V. The genetic signature of range expansion in a compilation of conserved polymerase chain reaction primers. Ann Entomol disease vector-the black-legged tick. J Hered. 2017;108:176–83. Soc Am. 1994;87:651–701. 42. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA 19. Ogden N. Changing geographic ranges of ticks and tick-borne pathogens: polymorphism data. Bioinformatics. 2009;25:1451–2. drivers, mechanisms and consequences for pathogen diversity. Front Cell 43. Excoffier L, Lischer HE. Arlequin suite ver 3.5: a new series of programs to Infect Microbiol. 2013;3:46. perform population genetics analyses under Linux and Windows. Mol Ecol 20. Criscione CD, Poulin R, Blouin MS. Molecular ecology of parasites: elucidating Resour. 2010;10:564–7. ecological and microevolutionary processes. Mol Ecol. 2005;14:2247–57. 21. Gooding RH. Genetic variation in arthropod vectors of disease-causing 44. Meirmans PG, Hedrick PW. Assessing population structure: FST and related organisms: obstacles and opportunities. Clin Microbiol Rev. 1996;9:301–20. measures. Mol Ecol Resour. 2011;11:5–18. 22. Gandon S, Michalakis Y. Local adaptation, evolutionary potential and host- 45. Rogers AR, Harpending H. Population growth makes waves in the parasite coevolution: interactions between migration, mutation, population distribution of pairwise genetic differences. Mol Biol Evol. 1992;9:552–69. size and generation time. J Evol Biol. 2002;15:451–62. 46. Harpending HC. Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Hum Biol. 1994;66:591–600. 23. Kubasu S. The ability of Rhipicephalus appendiculatus (Acarina: Ixodidae) 47. Fu YX. Statistical tests of neutrality of mutations against population growth, stocks in Kenya to become infected with Theileria parva. Bull Entomol Res. hitchhiking and background selection. Genetics. 1997;147:915–25. 1992;82:349–53. 48. Tajima F. Statistical method for testing the neutral mutation hypothesis by 24. Odongo DO, Ueti MW, Mwaura SN, Knowles DP, Bishop RP, Scoles GA. DNA polymorphism. Genetics. 1989;123:585–95. Quantification of Theileria parva in Rhipicephalus appendiculatus (Acari: Ixodidae) confirms differences in infection between selected tick strains. J 49. Kumar S, Stecher G, Tamura K. MEGA7: Molecular Evolutionary Genetics Med Entomol. 2009;46:888–94. Analysis Version 7.0 for bigger datasets. Mol Biol Evol. 2016;33:1870–4. 25. Young AS, Dolan TT, Mwakima FN, Ochanda H, Mwaura SN, Njihia GM, et al. 50. Leigh JW, Bryant D. popart: full-feature software for haplotype network Estimation of heritability of susceptibility to infection with Theileria parva in construction. Methods Ecol Evol. 2015;6:1110–6. the tick Rhipicephalus appendiculatus. Parasitology. 1995;111:31–8. 51. Speybroeck N, Madder M, Thulke HH, Mtambo J, Tirry L, Chaka G, et al. 26. Ochanda H, Young A, Medley G, Perry B. Vector competence of 7 Variation in body size in the tick complex Rhipicephalus appendiculatus/ rhipicephalid tick stocks in transmitting 2 Theileria parva parasite stocks Rhipicephalus zambeziensis. J Vector Ecol. 2004;29:347–54. from Kenya and Zimbabwe. Parasitology. 1998;116:539–45. 52. Mtambo J, Madder M, Van Bortel W, Chaka G, Berkvens D, Backeljau T. Further evidence for geographic differentiation in R. appendiculatus (Acari: 27. McCoy K. The population genetic structure of vectors and our Ixodidae) from Eastern and Southern provinces of Zambia. Exp Appl Acarol. understanding of disease epidemiology. Parasite. 2008;15:444–8. 2007;41:129–38. 28. Le Roux J, Wieczorek A. Molecular systematics and population genetics of biological invasions: towards a better understanding of invasive species 53. Excoffier L, Foll M, Petit RJ. Genetic consequences of range expansions. management. Ann Appl Biol. 2009;154:1–7. Annu Rev Ecol Evol Syst. 2009;40:481–501. 29. Mtambo J, Madder M, Van Bortel W, Geysen D, Berkvens D, Backeljau T. 54. Braverman JM, Hudson RR, Kaplan NL, Langley CH, Stephan W. The Genetic variation in Rhipicephalus appendiculatus (Acari: Ixodidae) from Zambia: hitchhiking effect on the site frequency spectrum of DNA polymorphisms. correlating genetic and ecological variation with Rhipicephalus appendiculatus Genetics. 1995;140:783–96. from eastern and southern Africa. J Vector Ecol. 2007;32:168–75. 55. Simonsen KL, Churchill GA, Aquadro CF. Properties of statistical tests of 30. Kanduma EG, Mwacharo JM, Mwaura S, Njuguna JN, Nzuki I, Kinyanjui PW, neutrality for DNA polymorphism data. Genetics. 1995;141:413–29. et al. Multi-locus genotyping reveals absence of genetic structure in field 56. Ray N, Currat M, Excoffier L. Intra-deme molecular diversity in spatially populations of the brown ear tick (Rhipicephalus appendiculatus) in Kenya. expanding populations. Mol Biol Evol. 2003;20:76–86. Ticks Tick Borne Dis. 2016;7:26–35. 57. Frankham R. Relationship of genetic variation to population size in wildlife. Conserv Biol. 1996;10:1500–8. 31. Kanduma EG, Mwacharo JM, Githaka NW, Kinyanjui PW, Njuguna JN, Kamau LM, et al. Analyses of mitochondrial genes reveal two sympatric but 58. Emerson BC, Paradis E, Thébaud C. Revealing the demographic histories of genetically divergent lineages of Rhipicephalus appendiculatus in Kenya. species using DNA sequences. Trends Ecol Evol. 2001;16:707–16. Parasit Vectors. 2016;9:353. 59. Robbertse L, Baron S, van der Merwe NA, Madder M, Stoltsz WH, Maritz- 32. Berkvens DL, Geysen DM, Chaka G, Madder M, Brandt JR. A survey of the Olivier C. Genetic diversity, acaricide resistance status and evolutionary ixodid ticks parasitising cattle in the Eastern Province of Zambia. Med Vet potential of a Rhipicephalus microplus population from a disease-controlled Entomol. 1998;12:234–40. cattle farming area in South Africa. Ticks Tick Borne Dis. 2016;7:595–603. Amzati et al. Parasites & Vectors (2018) 11:329 Page 19 of 19 60. Avise JC. Phylogeography: The history and formation of species. USA: Harvard University Press; 2000. 61. Slatkin M, Hudson RR. Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations. Genetics. 1991; 129:555–62. 62. Rogers AR, Fraley AE, Bamshad MJ, Watkins WS, Jorde LB. Mitochondrial mismatch analysis is insensitive to the mutational process. Mol Biol Evol. 1996;13:895–902. 63. Cristescu ME. Genetic reconstructions of invasion history. Mol Ecol. 2015;24:2212–25. 64. Biek R, Real LA. The landscape genetics of infectious disease emergence and spread. Mol Ecol. 2010;19:3515–31. 65. Ostfeld RS. Climate change and the distribution and intensity of infectious diseases. Ecology. 2009;90:903–5. 66. Manel S, Holderegger R. Ten years of landscape genetics. Trends Ecol Evol. 2013;28:614–21. 67. Berkvens DL, Pegram RG, Brandt JR. A study of the diapausing behaviour of Rhipicephalus appendiculatus and R. zambeziensis under quasi-natural conditions in Zambia. Med Vet Entomol. 1995;9:307–15. 68. Murrell A, Campbell NJH, Barker SC. Phylogenetic analyses of the rhipicephaline ticks indicate that the genus Rhipicephalus is paraphyletic. Mol Phylogenet Evol. 2000;16:1–7. 69. Burger TD, Shao R, Barker SC. Phylogenetic analysis of mitochondrial genome sequences indicates that the cattle tick, Rhipicephalus (Boophilus) microplus, contains a cryptic species. Mol Phylogenet Evol. 2014;76:241–53. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Parasites & Vectors Springer Journals

Mitochondrial phylogeography and population structure of the cattle tick Rhipicephalus appendiculatus in the African Great Lakes region

Free
19 pages

Loading next page...
 
/lp/springer_journal/mitochondrial-phylogeography-and-population-structure-of-the-cattle-vc6fTskjZX
Publisher
Springer Journals
Copyright
Copyright © 2018 by The Author(s).
Subject
Biomedicine; Parasitology; Entomology; Tropical Medicine; Infectious Diseases; Veterinary Medicine/Veterinary Science; Virology
eISSN
1756-3305
D.O.I.
10.1186/s13071-018-2904-7
Publisher site
See Article on Publisher Site

Abstract

Background: The ixodid tick Rhipicephalus appendiculatus is the main vector of Theileria parva, wich causes the highly fatal cattle disease East Coast fever (ECF) in sub-Saharan Africa. Rhipicephalus appendiculatus populations differ in their ecology, diapause behaviour and vector competence. Thus, their expansion in new areas may change the genetic structure and consequently affect the vector-pathogen system and disease outcomes. In this study we investigated the genetic distribution of R. appendiculatus across agro-ecological zones (AEZs) in the African Great Lakes region to better understand the epidemiology of ECF and elucidate R. appendiculatus evolutionary history and biogeographical colonization in Africa. Methods: Sequencing was performed on two mitochondrial genes (cox1 and 12S rRNA) of 218 ticks collected from cattle across six AEZs along an altitudinal gradient in the Democratic Republic of Congo, Rwanda, Burundi and Tanzania. Phylogenetic relationships between tick populations were determined and evolutionary population dynamics models were assessed by mismach distribution. Results: Population genetic analysis yielded 22 cox1 and 9 12S haplotypes in a total of 209 and 126 nucleotide sequences, respectively. Phylogenetic algorithms grouped these haplotypes for both genes into two major clades (lineages A and B). We observed significant genetic variation segregating the two lineages and low structure among populations with high degree of migration. The observed high gene flow indicates population admixture between AEZs. However, reduced number of migrants was observed between lowlands and highlands. Mismatch analysis detected a signature of rapid demographic and range expansion of lineage A. The star-like pattern of isolated and published haplotypes indicates that the two lineages evolve independently and have been subjected to expansion across Africa. (Continued on next page) * Correspondence: gastonamzati@gmail.com Unit of Integrated Veterinary Research, Department of Veterinary Medicine, Faculty of Sciences, Namur Research Institute for Life Sciences (NARILIS), University of Namur (UNamur), Rue de Bruxelles 61, 5000 Namur, Belgium Research Unit of Veterinary Epidemiology and Biostatistics, Faculty of Agricultural and Environmental Sciences, Université Evangélique en Afrique, P.O. Box 3323, Bukavu, Democratic Republic of the Congo Full list of author information is available at the end of the article © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Amzati et al. Parasites & Vectors (2018) 11:329 Page 2 of 19 (Continued from previous page) Conclusions: Two sympatric R. appendiculatus lineages occur in the Great Lakes region. Lineage A, the most diverse and ubiquitous, has experienced rapid population growth and range expansion in all AEZs probably through cattle movement, whereas lineage B, the less abundant, has probably established a founder population from recent colonization events and its occurrence decreases with altitude. These two lineages are sympatric in central and eastern Africa and allopatric in southern Africa. The observed colonization pattern may strongly affect the transmission system and may explain ECF endemic instability in the tick distribution fringes. Keywords: East Coast fever, Theileria parva, Rhipicephalus appendiculatus, Phylogenetic, cox1, 12S rRNA, Ticks, Evolutionary history, Agro-ecological zones, Population genetics Background The adaptive mechanism of R. appendiculatus to chan- The ixodid brown ear tick Rhipicephalus appendiculatus ged environment suggests genetic divergence between is the main vector of the protozoan pathogen Theileria geographical stocks and phenotypic variations, including parva, the causative agent of a fatal lymphoproliferative diapause behaviour and vector competence to transmit T. cattle disease known as East Coast fever (ECF). East Coast parva [23, 24]. The degree to which different R. appendi- fever is a highly pathogenic and the most economically culatus stocks acquire and transmit T. parva is partially important tick-borne disease of cattle in 12 sub-Saharan genetically dependent because of the heritability of their African countries, including Burundi, Democratic Repub- susceptibility to infection [25]. Furthermore, there is a lic of Congo and Rwanda [1]. Rhipicephalus appendicula- similar extensive genetic variation among T. parva strains tus is the most abundant tick in the Great Lakes region of in the field associated with different clinical features and Central Africa, where its burden and distribution vary sig- disease outcomes, and variable cross-immunity levels. nificantly among agro-ecological zones (AEZs) [2–4]. The Studies suggest that there are also specific interactions geographical dispersal pattern and population dynamics of between R. appendiculatus and T. parva stocks in the this tick are mainly driven by climatic conditions, vegeta- transmission dynamic system, significantly affecting the tion, host availability and mobility, grazing system and epidemiology of ECF [26]. Thus, phylogenetic and management practices [5, 6]. ecological analyses should provide useful information to The Great Lakes region of Central Africa is charac- control ECF by determining: (i) the genetic structure of R. terised by cattle movement within and between countries appendiculatus populations; (ii) its dispersal pattern; and for trade, breeding and pasture [7, 8]. During the (iii) its demographic history [27, 28]. pre-colonial period, immigrants originating from Rwanda The genetic diversity of R. appendiculatus has been and Burundi settled with their cattle in Congo in search of studied in different African countries using various mo- grazing lands. In addition, political unrest during the past lecular tools, such as mitochondrial DNA [14, 29] and three decades and rapidly increasing demand for animal micro- and minisatellite DNA markers [30]. Two distinct products increased significantly the importation of live an- genetic groups have been described, namely the eastern imals [9, 10]. This cross-border movement of cattle across and the southern African lineages [29]. More recently, geographical areas, together with bioclimatic conditions Kanduma et al. [31] found that the two genetic groups suitable for R. appendiculatus, could play a significant role are present in Kenya. These evidences show that R. in spreading ticks and pathogens [11–14]. Therefore, the appendiculatus genetic groups may have a wide geo- spread and establishment of ticks from one geographical graphical range, with different ecological preferences region to another might be setting up a complexity in the and phenology in sub-Saharan Africa [32–34], due to epidemiological status and control of the disease they differences in body size [35] and diapause induction and transmit [15–17]. Thus, predicting vector-borne pathogen intensity [36, 37]. dynamics and emergence relies on better understanding Major gaps in knowledge still exist concerning the of mechanisms underlying the genetic structure of their agro-ecological colonization and establishment of the R. vectors [18, 19]. appendiculatus lineages in the Great Lakes region of Ecological establishment and population genetic struc- Central Africa, where cattle mobility seems to be the main ture of ticks can be affected by founder events and gene factor of tick dispersal and epidemic instability of ECF [2, flow, largely due to their dispersal across geographical 38]. Thus, further studies on the population structure and zones through host migration [18, 20]. Arthropod vectors phylogeography of R. appendiculatus, including ticks from then respond differently to evolutionary forces such as mi- distinct agro-ecological conditions of DR Congo, Burundi gration, mutation, selection and genetic drift (promoted and Rwanda, are important to shed light on the intra and by bottlenecks) in their new environment [21, 22]. inter population variation, the dispersal pattern and the Amzati et al. Parasites & Vectors (2018) 11:329 Page 3 of 19 historical dynamics of the characterised lineages in various close to the equator, the wide range of altitudes (700 to ecological situations of Africa. The objective of this study 3000 m above sea level) mitigates significantly the was to assess the evolutionary relationship between R. tropical climate attributes and therefore the AEZs limits. appendiculatus populations and to investigate the impact Generally, the rainfall period is long and bimodal (with of geographical locations on its genetic structure, to better variations in the duration between AEZs). An early rainy understand the epidemiology of ECF in the Great Lakes season starts approximately in September and ends in region. To achieve this objective, we sequenced and ana- December, while the late rainy season lasts from February lysed fragments of the cytochrome c oxidase subunit 1 to May, followed by a dry season from June to August. (cox1) and the 12S ribosomal RNA (12S rRNA) gene loci. There is an intervening short dry period between the two The genetic structure of R. appendiculatus provides clues rainy seasons of about 15 days in January and/or February to a better understanding of the epidemiology of ECF and depending on the AEZs. Rainfall and temperature are insight into the development of targeted selective manage- strongly influenced by altitude ranges. Rainfall increases ment and control strategies. while temperature decreases with increasing elevation. The average temperature ranges between 12–25 °C and Methods the annual average rainfall from 800–1100 mm in the Study area Ruzizi Valley to 1300–2000 mm in the highlands. Ticks were collected from cattle in six agro-ecological In eastern DRC, ticks were collected along an altitudinal zones (AEZs) of the Central African Great Lakes region transect in South-Kivu province. Based on elevation and (Democratic Republic of Congo, Rwanda and Burundi) geographical position within the transect, three major AEZs (Fig. 1). The three countries share the Ruzizi Valley, were defined from lowlands to highlands, namely, DRC which consists of lowlands along the Ruzizi River flow- AEZ1 (lowlands) in the Ruzizi Valley, DRC AEZ2 (mid- ing between Lake Kivu and Lake Tanganyika. A sum- lands) in the administrative district of Walungu and DRC mary description of the characteristics of the six AEZs is AEZ3 (highlands) in the district of Kabare and the high- presented in Table 1. Briefly, although the study area is lands part of Walungu (Mulumemunene). The lowlands Fig. 1 Sampling sites of Rhipicephalus appendiculatus ticks in DRC, Rwanda and Burundi. a Map of Africa showing the study area (in grey) and other countries where the tick was previously sequenced (indicated by their names). b Sampling localities of R. appendiculatus and their altitudes (squares: AEZ1 altitude < 1200 m; circles: AEZ2 altitude 1200–1600 m; and triangles: AEZ3 altitude > 1600 m). The sites represented by empty circles and triangles show sampling locations described by Mtambo et al. [29] Amzati et al. Parasites & Vectors (2018) 11:329 Page 4 of 19 Table 1 Geographical and climatic attributes of the six agro-ecological zones (AEZ) Country Agro-ecological zone (AEZ) Altitude (m) Temperature (°C) Rainfall (mm/year) Rainy season Sample size (no. of ticks) DRC AEZ1 780–1100 23–25 800–1000 October-April 46 AEZ2 1200–1600 17–21 1000–1500 September-May 54 AEZ3 1600–2800 12–19 1350–2000 September-May 46 Burundi AEZ1 774–1100 23–24 800–1100 November-May 26 AEZ3 1700–2800 14–15 1300–2000 September-May 17 Rwanda AEZ2 1200–1500 21–24 800–950 November-May 20 Notes: DRC (AEZ1: Lowlands, AEZ2: Midlands, AEZ3: Highlands); Burundi (AEZ1: Lowlands, AEZ3: Highlands); Rwanda (AEZ2: eastern low plateau which is the lowlands of Rwanda as described by Bazarusanga et al. [3]) Sequences previously described in Mtambo at al. [29] are not included in the 20 samples from Rwanda and thus were not used in the population genetic analysis presented in Tables 2, 3 and 4 region (DRC AEZ1) is the warmest AEZ, characterised by a obtained from the Sokoine University of Agriculture in tropical dry climate (semi-arid), with a cool dry season of Tanzania. 4–5 months. Cattle are raised generally in an open grazing system in communal pastures of savannah grassland. In DNA extraction, PCR amplification and sequencing contrast, the “upland Kivu” (DRC AEZ2 and AEZ3) has a Total genomic DNA was isolated from whole individual montane humid tropical climate and is much cooler. The ticks using the DNeasy® Blood and Tissue Kit (Qiagen warm dry season lasts for 3–4 months, with occasional and GmbH, Hilden, Germany) according to the standard poorly distributed rainfall. manufacturer’s protocol, except that an additional incu- Ticks from Burundi were collected from two AEZs. bation of 10 min at 56 °C was added after mixing the The sampled administrative districts included Rugombo sample with 200 μl buffer AL to ensure optimal cell lysis. and Gihanga in the Imbo lowlands (Burundi AEZ1) and Prior to extract DNA, ticks were washed twice in double Muramvya, Mwaro and Mugamba districts of the distilled water and left to dry for 10 min at room Congo-Nile highlands (Burundi AEZ3). The Imbo region temperature and homogenized by grinding. extends along the Ruzizi River and the North of Lake Given the suitability of mitochondrial genes to dis- Tanganyika. The vegetation is composed of savannah criminate intraspecific variation in ticks [29, 31], we and small patches of forest. amplified the cox1and 12S rRNA gene loci to assess the Ticks from Rwanda were obtained from the eastern low genetic diversity and phylogenetic relationships of R. plateau in Nyagatare, Gatsibo and Kayonza districts appendiculatus populations. A 710 bp fragment of cox1 (Rwanda AEZ2). This region includes the savannah of gene locus was amplified with the forward primer eastern province of Rwanda. The eastern semi-arid plateau LCO1490 (5'-GGT CAA CAA ATC ATA AAG ATA zone is the most important pastoral region in Rwanda, TTG G-3') and the reverse primer HC02198 (5'-TAA holding 40% of the national herd raised in a communal ACT TCA GGG TGA CCA AAA AAT CA-3') as previ- grazing system. The vegetation type is largely savannah ously described by Folmer et al. [40]; for the 12S rRNA and some river bank woods. gene locus we used the primers SR-J-1499 (5'-TAC TAT GTT ACG ACT TAT-3') and SR-N-14594 (5'-AAA CTA Tick sampling and morphological identification GGA TTA GAT ACC C-3') with a fragment size of 380 Ticks were collected from cattle during a cross-sectional bp, described in Simon et al. [41]. PCR amplifications of survey conducted from February to April 2015 (late both cox1 and 12S rRNA gene fragments were per- rainy season) in the six AEZs. In each AEZ, 8 to 12 cat- formed using 50 ng of genomic DNA, 25 μl of 2× Accu- tle herds were selected randomly using a random num- Power® PCR PreMix (Bioneer PCR-PreMix, Seoul, South ber generator in Microsoft Excel. From ticks collected in Korea), 0.2 μM each of forward and reverse primers, and each herd representing a location or village, 4 to 5 ticks nuclease free water added up to a final reaction volume were selected using the same random process for further of 50 μl. The thermal cycling program for cox1 consisted analysis. The number of ticks sampled in each popula- of an initial denaturation at 95 °C for 3 min followed by tion is shown in Table 1. All ticks were directly 35 cycles of denaturation at 94 °C for 1 min, annealing immersed in 70% ethanol for preservation and morpho- at 40 °C for 1 min and extension at 72 °C for 1 min. The logically identified using the identification key of Walker final extension was carried out for 10 min at 72 °C. PCR et al. [39]. Morphological identification was confirmed at parameters of 12S rRNA gene fragment were as de- the tick unit at the International Livestock Research In- scribed for cox1, except that the annealing temperature stitute (ILRI, Kenya). Ten additional specimens originat- was 52 °C and the extension time was 90 s. PCR prod- ing from Simanjiro district in northern Tanzania were ucts were analysed by electrophoresis on a 1.8% agarose Amzati et al. Parasites & Vectors (2018) 11:329 Page 5 of 19 gel. Amplicons were purified using the QIAquick® PCR migrants. These analyses were performed for combined Purification Kit (Qiagen GmbH, Hilden, Germany) fol- data (of the main haplogroups identified) to understand lowing the manufacturer’s instructions. Both strands the effect of coexistence of the haplogroups in the gen- were sequenced directly with the same primers used for etic structure and diversity. PCR, on an ABI 3730 capillary sequencer (Applied Bio- systems, California, USA). Demographic and spatial expansion history analyses The historical population dynamics and structure was Sequences editing, blasting and multiple alignments inferred from mismatch distribution of cox1 haplotypes Forward and reverse chromatograms for each individual implemented in Arlequin [45], which compared the ob- tick were visually checked. Sequences were manually served distribution of pairwise nucleotide differences be- edited and assembled using CLC Main Workbench soft- tween haplotypes and that expected under a population ware v7.8.1 (CLC Bio, Aarhus, Denmark). Multiple se- expansion model for each population, haplogroup and quences were aligned with CrustalW algorithm using overall data. Multimodal mismatch pattern is assumed default settings in the same software. The sequences to define a population of demographic equilibrium or were then trimmed to exclude poor quality bases and constant size, whereas sudden or expanding population obtain uniform sizes. The final sequence sizes were 586 is characterised by unimodal distribution. The sum of bp and 346 bp for cox1 and 12S rRNA genes, respect- square deviation (SSD) were estimated to determine if ively. Aligned and trimmed sequences were subjected to the observed mismatch deviated significantly from a a BLAST search against all NCBI nucleotide databases, population expansion model, and the Harpending’s rag- with default settings to confirm their species identity. A gedness index (RI) were used to evaluate the smoothness total of 219 sequences for cox1 and 126 for 12S were ob- of mismatch distribution [46]. In addition, to detect de- tained after quality check processing. To check for amp- viation from neutrality expectations or mutation-drift lification of nuclear pseudogenes and confirm protein equilibrium we performed analysis of Fu’s Fs [47] and coding, all cox1 nucleotide sequences were translated Tajima’s D [48] statistics in Arlequin and DnaSP. Taji- into their amino acid sequences to examine the presence ma’s D test is based on the difference between the num- of ambiguous stop codons for correct coding of inverte- ber of polymorphic sites and the mean number of brate mitochondrial DNA. nucleotide pairwise differences, while Fu’s Fs test is based on haplotype frequencies. The significance of all Genetic diversity and population genetic structure these statistics was tested by bootstrap resampling of Multiple sequences extracted from the alignments were 1000 coalescent simulations. Significant negative values used to construct haplotypes in DnaSP software v5.10.01 of neutrality statistics should indicate a signature of his- [42]. Genetic variation within populations was estimated torical event of population expansion, whereas signifi- on cox1 gene sequences. Population genetic indices rep- cantly positive values indicate events such as recent resented by number of haplotypes (h), number of segre- population bottleneck, population subdivision or pres- gating sites (S), haplotype diversity (Hd), mean number ence of some recent immigrants in a population. Values of pairwise nucleotide differences within population (K) near zero and not significant, indicate that population and nucleotide diversity (π) were calculated for each size is constant. AEZs (named populations) and for the overall data set using Arlequin v3.5.2.2 [43] and DnaSP. The population Phylogenetic and phylogeographical reconstruction genetic structure among AEZs and among haplogroups Different published haplotype sequences of R. appendi- was evaluated by an analysis of molecular variance culatus from South Africa, Kenya, Grande Comore, (AMOVA) performed in Arlequin. The significance of Zambia, Zimbabwe, Uganda and Rwanda were retrieved AMOVA fixation indices was evaluated based on 1,023 from the National Center for Biotechnology Information random permutations. In addition, to assess the level of (NCBI) database (see Additional file 1: Table S1) and genetic distance/differentiation between populations, we were compared with sequences obtained in the present estimated gene flow expressed as of the absolute number study. Phylogenetic reconstruction was performed separ- of migrants (Nm) exchanged among populations, average ately on cox1 and 12S rRNA gene sequences to deter- number of nucleotide differences between populations mine the relationship among populations and possible (K ) and population pairwise genetic differentiation historical dispersal event. To find the evolutionary sub- XY (F ) also in Arlequin [44]. Their significance was tested stitution model that best describe the evolution of cox1 ST using 1000 random permutations. The genetic differenti- and 12S rRNA gene sequences, we performed a hier- ation and population structure statistics were tested archical likelihood ratio test, based on the lowest Bayes- under the hypothesis that different populations are rep- ian information criterion using MEGA v7.0 [49]. The resented by distinct genetic groups or are exchanging nucleotide substitution model selected was then used to Amzati et al. Parasites & Vectors (2018) 11:329 Page 6 of 19 perform a Neighbor-Joining (NJ) and/or Maximum Like- molecular identity (99–100%) with known sequences of lihood (ML) algorithm in MEGA. The stability and R. appendiculatus found in GenBank (Additional file 2: branches support were obtained using 1000 bootstrap Table S2). Haplotype sequences (55 and 23 sequences permutations. Rhipicephalus eversti and Rhipicephalus for cox1and 12S rRNA, respectively) obtained for microplus from this study (GenBank accession numbers each of the six AEZs were deposited and are avail- MF458972 and MF458973 for cox1 and MF479198 and able in GenBank (GenBank accession numbers: MF479199 for 12S RNA genes) and Rhipicephalus tura- MF458895-MF458949 and MF479166-MF479188 for nicus obtained from the GenBank (accession numbers cox1and 12S rRNA genes, respectively). KU880574 and DQ849231 for cox1and 12S rRNA genes, respectively) were used as outgroup taxa. A Phylogenetic relationships and haplotypes distribution Median-Joining (MJ) network was constructed to inves- The overall sequence analysis revealed that cox1had 27 tigate the phylogenetic and ancestral relationship among polymorphic sites, 21 of which were parsimony inform- haplotypes using PopArt Software [50]. ative and 6 were singletons, yielding 22 haplotypes (Additional file 3:Table S3). cox1 amino acid sequences Results did not contain any internal stop codon or indel. Most Morphological and molecular identification of ticks nucleotide mutations were synonymous, except one PCR fragments of the mitochondrial cox1 and 12S rRNA non-synonymous change identified at position 32 of the gene loci were successfully generated from 219 and 126 haplotype CH22 (change from an Alanine to a Threonine). individual ticks, respectively, originating from the three The highest number of cox1 haplotypes was found in AEZs of DRC, the two AEZs of Burundi, one AEZ of DRC AEZ1, while the lowest was observed in Burundi Rwanda and specimens from Tanzania. The 10 se- AEZ3 (Table 2). The 12S rRNA gene provided 10 poly- quences from Tanzania and the sequences previously de- morphic sites, five of which were parsimony informative, scribed in Rwanda were not included in the population defining 9 haplotypes. The 22 cox1 haplotype sequences genetic and structure analyses. Generated nucleotide se- obtained in this study were submitted to GenBank under quences were aligned with the reference haplotype se- accession numbers MF458950-MF458971 and the nine quences retrieved from the GenBank (Additional file 1: 12S rRNA gene haplotypes were deposited under acces- Table S1). The final fragment length obtained for cox1 sion numbers MF479189-MF479197. The phylogenetic re- was 586 bp and 12S rRNA yielded a fragment of 346 bp lationships among cox1haplotypes inferredbya NJ long, with no indels detected in both genes. Their mor- phylogenetic tree (Fig. 2), a ML tree (Fig. 3) and a MJ net- phological identification has been confirmed by the high work (Fig. 4) produced identical topologies and detected Table 2 Rhipicephalus appendiculatus cox1 haplotypes distribution (%) and genetic diversity indices in six agro-ecological zones of the Democratic Republic of Congo, Burundi and Rwanda a,b c Country AEZ n h Haplotype (frequency in %) Haplogroup Genetic diversity indices (%) AB S Hd (SD) K π (SD) (PIS) DRC AEZ1 46 10 CH1(30), CH2(28), CH5(6), CH6(2), CH7(4), CH8(2), CH11(9), CH13(11), 85 15 17 0.82 4.4 (2.2) 0.007 CH16(4), CH17(2) (16) (0.03) (0.002) AEZ2 54 12 CH1(18), CH2(30), CH5(28), CH7(2), CH11(2), CH12(4), CH13(4), 92 8 20 0.81 3.1 (1.6) 0.005 CH18(4), CH19(4), CH20(2), CH21(2), CH22(2) (17) (0.03) (0.001) AEZ3 46 7 CH1(17), CH2(33), CH5(28), CH11(9), CH12(6), CH13(4), CH14(2) 96 4 16 0.79 2.4 (1.3) 0.004 (15) (0.03) (0.001) Burundi AEZ1 26 8 CH1(50), CH2(19), CH3(4), CH4(4), CH5(11), CH6(4), CH7(4), CH10(4) 96 4 18 0.72 2.01 (1.2) 0.003 (3) (0.08) (0.001) AEZ3 17 5 CH1(29), CH2(35), CH5(23), CH8(6), CH9(6) 100 0 4 (2) 0.77 1.1 (0.9) 0.002 (0.06) (0.0003) Rwanda AEZ2 20 8 CH1(30), CH2(20), CH4(5), CH5(5), CH7(20), CH11(5), CH13(10), 70 30 18 0.85 6.3 (3.1) 0.011 CH15(5) (13) (0.05) (0.002) Total 209 22 – 90 10 27 0.81(0.01) 3.4 (1.7) 0.006 (21) (0.0006) Abbreviations: AEZ agro-ecological zones, n number of sequences, h number of haplotypes, S segregation sites, PIS parsimony informative sites, Hd haplotype diversity, SD standard deviation; K mean number of pairwise nucleotide differences, π nucleotide diversity, CH1-22 names of cox1 haplotypes Haplotypes belonging to the haplogroup B are underlined Bold indicates shared haplotypes by all agro-ecological zones A and B are haplogroup names Amzati et al. Parasites & Vectors (2018) 11:329 Page 7 of 19 Fig. 2 Phylogenetic tree of R. appendiculatus cox1 haplotypes. The evolutionary history was inferred by using the neighbor-joining method based on the Tamura 3-parameter model. A discrete Gamma distribution was used to model evolutionary rate differences among sites. Bootstrap values (> 80) are displayed above nodes. CH1-22 are haplotype names. The values in parentheses correspond to the frequency of each haplotype. KU725893 and AF132833 are GenBank accession numbers for R. appendiculatus sequences used as reference haplotypes from Kenya and Zimbabwe, respectively. Two species (R. eversti and R. microplus) obtained in this study and R. turanicus from GenBank (accession number: KU880574) were included as outgroups two distinct clades or lineages strongly-supported by a NJ The distribution of haplotypes presented in Table 2 bootstrap value of 100%. The two lineages diverged at showed that there were shared and private cox1 haplo- least by 12 mutational steps (Additional file 3: Table S3) types confined to restricted AEZs. Three haplotypes but shared a wide range of agro-ecological and geograph- CH1 (27%), CH2 (28%) and CH5 (19%) were shared be- ical conditions in the Great Lakes region. The first lineage, tween all AEZs and were by far the most ubiquitous in named “haplogroup A”, was represented by the most fre- the region, accounting for 74% (154 sequences) of the quent haplotypes (CH1, CH2 and CH5) and comprised in overall dataset (Additional file 3: Table S3). Haplotype total 19 haplotypes consisting of 189 of the 209 sequences CH1 was detected in 13 (50%) of the 26 sequences from analysed (90%), whereas the second lineage “haplogroup lowlands of Burundi (Burundi AEZ1). Two haplotypes B” included three haplotypes (CH7, CH13, CH20) and (CH11 and CH13) defined by 10 and 11 sequences, re- had a total of 20 sequences (10%) (Fig. 2, Additional file 3: spectively, were common in all the AEZs of DRC and in Table S3). Haplogroup B comprised of haplotypes present Rwanda AEZ2. Haplotype CH12 was exclusive to DRC in most AEZs except the highlands region of Burundi AEZ2 and AEZ3. Nine out of the 22 haplotypes were (Burundi AEZ3). These cox1 haplogroups presented a found in single individuals; therefore, they belonged to star-like pattern in the MJ network with less frequent and particular AEZs (Additional file 3: Table S3). Further- single haplotypes connected together to the predominant more, two 12S rRNA haplotypes (12SH1 and 12SH2) or ancestral haplotypes generally by single substitutions were the most abundant, representing 90% of the 126 (Fig. 4), supporting the hypothesis of a recent population sequences analysed. Haplotypes 12SH4 and 12SH5 had expansion. The phylogenetic relationships found in the together 8 (6%) out of the 126 analysed sequences. The cox1 gene were fully congruent with those revealed by the presence of single haplotypes indicates high frequency of ML tree performed on 12S rRNA haplotypes. rare alleles, which suggests a recent population expansion. Amzati et al. Parasites & Vectors (2018) 11:329 Page 8 of 19 Fig. 3 Phylogenetic tree of cox1 haplotypes displaying the relationship between the R. appendiculatus specimens in sub-Saharan African countries. The evolutionary history was inferred by using the maximum likelihood method based on the Tamura 3-parameter model. A discrete Gamma distribution was used to model evolutionary rate differences among sites [5 categories (+G, parameter = 0.39)]. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. Bootstrap scores > 70 are displayed to support nodes. The values in bracket behind haplotype names correspond to the frequency of each haplotype. Haplotype sequences (CH1-22) obtained in the present study are indicated by a black square. Rhipicephalus eversti and R. microplus obtained in this study and R. turanicus (GenBank number: KU880574) were used as outgroups Population genetic diversity from 0.72 in Burundi AEZ1 to 0.85 in Rwanda AEZ2 and Population genetic indices were estimated using cox1nu- the nucleotide diversity (π) values ranged from 0.002 in cleotide sequences and are shown in Table 2. The overall Burundi AEZ3 to 0.011 in Rwanda AEZ2. The average haplotype diversity (Hd) and nucleotide diversity (π)were number of nucleotide pairwise differences (k) was 3.4 for 0.81 and 0.006, respectively. The haplotype diversity ranged the overall dataset, with the highest value observed in Amzati et al. Parasites & Vectors (2018) 11:329 Page 9 of 19 Fig. 4 Median-joining network of the 36 cox1 haplotypes for R. appendiculatus ticks across sub-Saharan African countries. Lines represent mutations and the dot corresponds to a possible intermediate haplotype. Each circle denotes a unique haplotype. Haplotype frequencies are not shown here, but their occurrences across Africa are presented in Table 5 Rwanda AEZ2 (6.3). Altogether, these population genetic (K ), the pairwise genetic distance (F ), and the number XY ST indices showed that the diversity of R. appendiculatus is of migrants (Nm) showed evidence of differentiation quite different among AEZs. Ticks from Rwanda AEZ2 between some R. appendiculatus populations. The popu- were more highly diverse than those from three AEZs of lation pairwise genetic distance (F ) varied from a nega- ST DRC and the two AEZs of Burundi. We also observed an tive and not significant value (F = -0.014, P = 0.87) with ST excess of singleton mutations in Burundi AEZ1 (15 out of infinite value of Nm (between DRC AEZ2 and AEZ3) to 18 polymorphic sites), suggesting a suddenpopulationex- the strongest differentiation (F =0.19, P =0.005)with ST pansion. These data analysed separately by haplogroups low number of migrants (Nm = 2.1) between Burundi (Table 3), showed that the differences of genetic diversity AEZ3 and Rwanda AEZ2. Pairwise F values were found ST among AEZs was largely affected by the frequency distribu- to be low and not significant between DRC AEZ3 and tion of the two cox1 haplogroups (Table 2). Burundi AEZ3 with infinite Nm and between DRC AEZ2 and Burundi AEZ3 associated to a high number of mi- Population structure and ecological differentiation based grants (Nm = 11,875), indicating an excess of gene flow on cox1 haplotypes between these zones. DRC AEZ1 did not differ significantly Analysis of molecular variance (AMOVA) based on cox1 with Burundi AEZ1 (F =0.022, P =0.19) andRwanda ST sequences revealed statistically significant variance among AEZ2 (F = 0.018, P = 0.24). Tick populations from ST the six AEZs analysed together for both combined Rwanda AEZ2 showed a strong genetic differentiation with haplogroups and haplogroup A alone (6%, P <0.001) those from DRC AEZ3 (F =0.18, P <0.001)and DRC ST (Additional file 4: Table S4). The genetic variability among AEZ2 (F =0.14, P = 0.03), indicating reduced or lack of ST individuals within AEZs explained the large majority of gene flow between these populations. These results were molecular variance (94%) for the overall dataset, suggest- confirmed by inter-population nucleotide differences (K ), XY ing an admixture between different populations and the which was significant when populations were significantly coexistence of the two genetically divergent lineages (Fig. 2, differentiated by the pairwise F statistic. Tick populations ST Table 2). This is consistent with moderate genetic struc- were compared by analysing haplogroup A sequences ture of R. appendiculatus in the Great Lakes region. The alone. The findings showed that ticks from Rwanda AEZ2 population differentiation indices were then calculated to and the two AEZs of Burundi were not genetically different. compare the genetic variation observed among AEZs They shared more migrants belonging to haplogroup A (Table 4). Both differentiation statistics based on pairwise (Table 4). In general, the population differentiation was ob- estimates of the Inter-population nucleotide differences served between lowlands and highlands AEZs. Amzati et al. Parasites & Vectors (2018) 11:329 Page 10 of 19 Table 3 Genetic diversity and evolutionary dynamics of the two haplogroups (A and B) identified from cox1 sequences of R. appendiculatus Genetic indices and Haplogroup A Haplogroup Overall data statistics B set DRC Burundi Rwanda Haplogroup AEZ2 A (overall) AEZ1 AEZ2 AEZ3 AEZ1 AEZ3 Diversity indices Number of sequences 39 50 44 25 17 14 189 20 209 Number of polymorphic 79 68 47 18 2 27 sites Number of haplotypes 8 9 6 7 5 6 19 3 22 Haplotype diversity (SD) 0.76 (0.044) 0.78 (0.032) 0.77 (0.033) 0.7 (0.084) 0.77 (0.057) 0.77 (0.089) 0.77 (0.016) 0.56 (0.063) 0.81 (0.01) Nucleotide diversity (SD) 0.002 0.0025 0.002 0.002 0.002 0.002 0.002 0.001 0.006 (0.0003) (0.0003) (0.0003) (0.0004) (0.0003) (0.0006) (0.0001) (0.0002) (0.0006) Neutrality tests Tajima's D (P-value) -0.6 (0.35) -0.79 (0.23) 0.096 (0.6) -1.4 -0.21 (0.45) -1.4 (0.07) -1.5 (0.032)* 0.24 (0.67) -0.93 (0.23) (0.049)* Fu's Fs (P-value) -2.4 (0.08) -2.5 (0.097) -0.1 (0.49) -2.6 -1.1 (0.18) -2.1 -10.4 0.2 (0.49) -3.8 (0.15) (0.027)* (0.041)* (0.001)* Demographic expansion Sum of Squared deviation 0.002 (0.66) 0.004 (0.31) 0.002 (0.49) 0.002 (0.72) 0.02 (0.17) 0.001 (0.9) 0.0008 0.029 0.017 (0.1) (P-value) (0.44) (0.049)* Harpending’s Raggedness 0.06 (0.39) 0.079 (0.21) 0.06 (0.38) 0.05 (0.71) 0.16 (0.12) 0.06 (0.73) 0.057 (0.2) 0.21 (0.1) 0.049 (0.51) index (P-value) Spatial expansion Sum of Squared deviation 0.002 (0.5) 0.004 (0.21) 0.002 (0.4) 0.002 (0.67) 0.02 (0.1) 0.001 (0.9) 0.0008 0.029 0.036 (0.28) (P-value) (0.26) (0.011)* Harpending’s Raggedness 0.064 (0.39) 0.079 (0.21) 0.062 (0.43) 0.05 (0.75) 0.16 (0.12) 0.06 (0.72) 0.057 (0.22) 0.21 (0.12) 0.049 (0.74) index (P-value) Abbreviations: D Tajima’s neutrality statistic, Fs Fu’s neutrality statistic *Values are statistically significant at P < 0.05; significance was determined using 1000 coalescent simulations Demographic and dispersal dynamics of R. appendiculatus haplogroup A, we detected significant evidence of demo- Demographic and spatial dynamics inferred from pairwise graphic and spatial expansion events from the unimodal nucleotide differences revealed bimodal pattern for the mismatch distribution, together with significantly nega- total dataset (Fig. 5a). Tajima’s D and Fu’s Fs were negative tive values for Tajima’s D (D = -1.5, P = 0.032) and Fu’s but not significant (Table 3). However, the sum of square Fs (Fs = -10.4, P = 0.001) statistics. A good fit of sudden deviation (SSD) and raggedness index (RI) for both demo- population expansion was also observed in this hap- graphic and range expansion did not deviate significantly logroup, based on sum of squared deviation values that from that expected under expansion model. The negative were not significant in all the cases: demographic (SSD = values of neutrality tests and the non-significance of SSD 0.0008, P = 0.44) and range (SSD = 0.0008, P = 0.26) ex- and RI indices, suggest a sudden expansion of R. appendi- pansion, with no significant variation of the raggedness culatus populations in the Great Lakes region. The popu- index for both models (Table 3). In contrast, haplogroup lation dynamics history was also analysed separately for B did not display any significant signature of expansion each of the six AEZs (Fig. 6, Additional file 5: Table S5). from the selective neutrality tests. In addition, the ob- Most AEZs showed a bimodal mismatch distribution, served mismatch pattern for this haplogroup deviated except in Burundi AEZ3 where the mismatch pattern was significantly from that expected under population expan- unimodal (Fig. 6e). The observed bimodal pattern sug- sion scenario (SSD = 0.029, P = 0.049 and P = 0.011 for gested population subdivision as shown by the existence demographic and spatial expansion, respectively), imply- of two well resolved haplogroups (Fig. 2). ing that haplogroup B did not experience any historical The population dynamics were then inferred for each population expansion. This group is characterised either haplogroup separately (Table 3). A unimodal distribution by a constant population size (demographic equilibrium) was observed in both haplogroups (Fig. 5b, c). For the or had experienced a weak signal of population Amzati et al. Parasites & Vectors (2018) 11:329 Page 11 of 19 Table 4 Population genetic statistics for pairwise comparison of different populations of R. appendiculatus from sequences of cox1 gene. Values in parentheses represent the P-value statistics Population 1 Population 2 Haplogroup A and B Haplogroup A K F Nm K F Nm XY ST XY ST DRC AEZ1 DRC AEZ2 3.8 (0.044)* 0.036 (0.044)* 13.2 1.5 (< 0.001)* 0.11 (< 0.001)* 4 DRC AEZ1 DRC AEZ3 3.6 (0.032)* 0.057 (0.032)* 8.3 1.5 (0.001)* 0.1 (0.001)* 4.3 DRC AEZ1 Burundi AEZ1 3.3 (0.15) 0.022 (0.19) 22.2 1.2 (0.31) 0.002 (0.33) 243.6 DRC AEZ1 Burundi AEZ3 3 (0.059) 0.056 (0.087) 8.4 1.3 (0.087) 0.043 (0.082) 11.2 DRC AEZ1 Rwanda AEZ2 5.4 (0.3) 0.018 (0.24) 26.8 1.3 (0.78) -0.026 (0.79) ∞ DRC AEZ2 DRC AEZ3 2.7 (0.85) -0.014 (0.87) ∞ 1.4 (0.93) -0.016 (0.93) ∞ DRC AEZ2 Burundi AEZ1 2.7 (0.003)* 0.060 (0.011)* 7.9 1.5(< 0.001)* 0.13 (< 0.001)* 3.2 DRC AEZ2 Burundi AEZ3 2.1 (0.23) 0.00004 (0.43) 11876 1.3 (0.46) -0.012 (0.54) ∞ DRC AEZ2 Rwanda AEZ2 5.2 (0.003)* 0.14 (0.003)* 3.1 1.6 (0.006)* 0.11 (0.004)* 3.8 DRC AEZ3 Burundi AEZ1 2.4 (< 0.001)* 0.076 (0.006)* 6.1 1.5(< 0.001)* 0.14 (< 0.001)* 3.2 DRC AEZ3 Burundi AEZ3 1.8 (0.32) -0.005 (0.43) ∞ 1.3 (0.39) -0.009 (0.46) ∞ DRC AEZ3 Rwanda AEZ2 5.1 (< 0.001)* 0.18 (< 0.001)* 2.2 1.6 (0.009)* 0.11 (0.011)* 3.9 Burundi AEZ1 Burundi AEZ3 1.6 (0.074) 0.034 (0.08) 12.2 1.2 (0.077) 0.063 (0.073) 7.4 Burundi AEZ1 Rwanda AEZ2 4.8 (0.030)* 0.14 (0.03)* 3.1 1.2 (0.97) -0.039 (0.98) ∞ Burundi AEZ3 Rwanda AEZ2 4.6 (0.005)* 0.19 (0.005)* 2.1 1.3 (0.16) 0.044 (0.14) 10.9 Abbreviations: K average number of nucleotide differences between populations, F pairwise genetic distance F-statistic based on nucleotide sequences (Right’s XY ST fixation index), Nm number of migrants between populations *Values are statistically significant at P < 0.05; significance was determined using 1000 coalescent simulations. bottleneck that reduced its diversity. When analysing the each of the six AEZs and for Tanzania (Table 5). demographic dynamics for samples belonging to haplogroup Twenty-eight haplotypes had been already described in dif- A in each AEZ (Table 3), population expansion signal was ferent African countries, and eight haplotypes (CH3, CH9, confirmed in all AEZs by mismatch analyses exhibiting uni- CH10, CH14, CH16, CH18, CH20 and CH33) were newly modal distribution (Additional file 6:FigureS1).For thesix described in the present study. Most of these new haplo- AEZs, none of the statistical comparisons between the ob- types were less abundant (Additional file 3:Table S3)and served and the expected distributions rejected the sudden diverged from the common ancestral haplotypes generally and range expansion assumptions based on the raggedness by only one substitution (Fig. 4). Haplotypes CH1, CH7, index and the sum of squared deviation. The neutrality indi- CH11 were the most ubiquitous and shared wide geo- ces were generally negative but not significant, except in graphical distribution in affected African countries. CH1 Burundi AEZ1 where both Tajima’s D (D =-1.4, P = 0.049) haplotype was shared by Kenya, Eastern Zambia, DRC and Fu’s Fs (Fs =-2.6, P = 0.027) statistics showed significant (AEZ1, 2 and 3), Burundi (AEZ1 and 3), Tanzania and negative values and in Rwanda AEZ2 where Fu’s Fs was Rwanda (AEZ2 and GenBank sequences), while haplotype negative (Fs = -2.1) and significant (P = 0.041). According to CH7 was reported in Kenya, South Africa, Zimbabwe, the mismatch distribution and negative values for neutrality Grande Comore, Eastern and southern provinces of tests, the hypothesis of population growth and spatial expan- Zambia, DRC (AEZ1, 2 and 3), Burundi (AEZ1) and sion models could not be rejected in the six AEZs, which Rwanda (AEZ2). Haplotype CH11 was present in Kenya, was consistent with a model of sudden expansion for each Rwanda (AEZ2 and GenBank sequences), Comoros and population subdivision. DRC (AEZ1, 2 and 3). Eighteen haplotypes mostly with unique sequences were found to be restricted to Kenya Phylogeographical structure and have not been reported in any other country. In the To study the phylogeographical structure of R. appendicu- same way, haplotype CH23 was particular to Uganda. This latus in Africa, the haplotype sequences found in this study country did not share any haplotype with other countries. along with those retrieved from GenBank (Additional file 1: The ML phylogenetic tree reconstructed using the 105 se- Table S1) were used to reconstruct the phylogenetic using quences indicated that globally, the African tick R. appen- ML tree and the MJ network methods based on cox1and diculatus is consistently clustered into two groups 12S rRNA genes. Thirty-six cox1 haplotypes were identified (haplogroups A and B) well-supported by a ML bootstrap in 105 sequences including 50 haplotypes from GenBank value of 100% (Fig. 3). Our 19 haplotypes that had been and the 55 haplotype sequences obtained in our study for described for haplogroup A (Table 2,Additional file 3: Amzati et al. Parasites & Vectors (2018) 11:329 Page 12 of 19 Fig. 5 cox1 mismatch distribution pattern for two haplogroups of R. appendiculatus. a shows the overall mismatch distribution pattern for all AEZs, b and c show the mismatch distribution of nucleotide sequences in haplogroups A and B, respectively. The x-axis shows the number of pairwise differences between pairs of haplotype sequences and the y-axis shows their frequencies (in %). The observed frequencies are represented by solid histograms and the simulated mismatch distributions expected under demographic expansion (solid black line) and under spatial expansion (dotted black line). Simulated curves under range and demographic expansion have same pattern in these figures, they overlapped Table S3) formed one clade with 19 haplotypes from Kenya (19/29), all the 3 haplotypes from Rwanda (3/3), six haplotypes from Eastern province of Zambia (6/7), all the five haplotypes from Tanzania, whereas our three haplo- types of haplogroup B were clustered with 10 haplotypes from Kenya, the single haplotype from Grande Comore, all haplotypes from South Africa (3/3), Zimbabwe (3/3), Uganda (2/2), Southern province of Zambia (2/2), and one haplotype from Eastern province of Zambia. In addition to ML tree, the MJ network also revealed that R. appendi- culatus is divided into two main groups in Africa, sepa- rated by seven mutational steps (Fig. 4). Similar findings were confirmed by 12S rRNA gene. Our 23 12S rRNA individual haplotypes from each of the six AEZs were analysed together with the 29 sequences ob- tained from GenBank (Additional file 1: Table S1). Fourteen haplotypes were observed, two most common (12SH1 and 12SH5) and 12 minors (defined mostly by one sequence or restricted to particular country) (Additional file 7:Table S6). Haplotype 12SH1 was common in DRC, Burundi, Rwanda, Kenya and Eastern province of Zambia, while haplotype 12SH5 was present in DRC, Rwanda, Zimbabwe, Comoros, South Africa, Eastern Zambia and Kenya. Six new haplo- types were not found outside the Great lakes region (12SH3, 12SH4 and 12SH6-H9). The NJ phylogenetic resolved these 12S rRNA haplotypes into two haplogroups (haplogroup A and B) supported by 97% bootstrap value (Additional file 8: Figure S2). Their pattern was identical to that observed from cox1 haplogroups. However, the Ugandan haplotype (12SH10: GenBank AF150028) was clustered in haplogroup A, showing that the haplogroup A is also present in Uganda. Unfortunately, we did not find its corresponding cox1 sequence. Discussion This study analysed the intraspecific variation of mitochon- drial DNA to better understand how factors such as agro-ecological zones and anthropogenic movements of cat- tle may affect population genetic structure and population expansion history of the tick R. appendiculatus,the main vector of T. parva in sub-Saharan Africa. We expected evi- dence of R. appendiculatus population expansion, gene flow Amzati et al. Parasites & Vectors (2018) 11:329 Page 13 of 19 Fig. 6 cox1 mismatch distribution pattern for six populations of R. appendiculatus. a, b and c show the mismatch distribution pattern for R. appendiculatus from DRC AEZs (AEZ1, 2 and 3, respectively); d and e represent the mismatch pattern of ticks from Burundi AEZ1 and AEZ3, respectively; f depicts the mismatch distribution of ticks from Rwanda AEZ2. The x-axis shows the number of pairwise differences between pairs of haplotype sequences and the y-axis shows their frequencies in %. The observed frequencies are represented by solid histograms. Black full line represents the expected distribution under sudden expansion model, and dotted line represents the distribution simulated under spatial expansion model. Simulated curves under spatial and demographic expansion have same pattern in (d), and they overlapped and different colonization patterns of tick lineages, facilitated groups [29], corresponding to our haplogroups A and B, by the reported cattle mobility in the Great Lakes region. respectively. This genetic grouping fitted well with the phenotypical, physiological and ecological variation stud- Two R. appendiculatus lineages that are more variable in ies which have previously distinguished two major sub- lowlands than highlands occur in the Great Lakes region populations of R. appendiculatus in Africa [32–37, 51]. The 22 haplotypes identified by DNA polymorphic ana- These phenotypic and physiological variations are largely lysis of cox1 gene locus were clustered into two associated to agro-ecological and geographical subdivi- well-defined major groups, named haplogroup A (the sions. Tropical areas with prolonged and marked dry most frequent) and haplogroup B. Similar grouping were seasons are more suitable for larger sized ticks express- obtained with 12S rRNA analyses. The two haplogroups ing high intensity of diapause and displaying univoltine identified in the present study have been previously de- phenology, corresponding to “southern African group” scribed as “east African” and “southern African” genetic or haplogroup B. Equatorial areas with bimodal or Amzati et al. Parasites & Vectors (2018) 11:329 Page 14 of 19 Table 5 Rhipicephalus appendiculatus cox1 haplotypes and their distribution among agro-ecological zones of the Great lakes region and other sub-Saharan African countries Haplotype Haplotypes from GenBank: Country (original haplotype name and Present study Haplogroup GenBank number) a e CH1 Kenya (H2: KU725891, H4: KU725893, H6: KU725895) , Rwanda Burundi (AEZ1, AEZ3), DRC (AEZ1, AEZ2, AEZ3), A f e g (H3:DQ901360) Zambia-east (H4: KX276942 , H5: DQ859266 , Tanzania (TZ18, TZ10, TZ08, TZ20), Rwanda f g H3:DQ901361 , H2: DQ859265 ) (AEZ2) CH3 – Burundi AEZ1 A CH4 Rwanda (H6: DQ901362) Burundi AEZ1, Rwanda (AEZ2) A b e CH6 Kenya (H5: KU725894) Burundi (AEZ1, AEZ3), DRC (AEZ1, AEZ2) A c e e e CH7 Kenya (H1: KU725890) , South Africa (H1: KX276939 , H1: KX276940 , H1: Burundi AEZ1, DRC (AEZ1, AEZ2, AEZ3), Rwanda B f g e DQ901356 ), Zambia-east (H1: DQ859261) , Zambia-south (H1:KX276943 , (AEZ2) g h i e H1: DQ859262) , Zimbabwe (AF132833 , KC503257 , H1: KX276944 ), Grande Comore (H1) CH9 – Burundi AEZ3 A CH10 – Burundi AEZ1 A d e f j CH11 Kenya (H3: KU725892) , Rwanda (H3: DQ901363) , Grande Comore (H3) DRC (AEZ1, AEZ2, AEZ3), Rwanda (AEZ2) A CH12 Kenya (H11: KU725900) DRC (AEZ2, AEZ3) A CH14 – DRC AEZ3 A CH16 – DRC AEZ1 A CH18 – DRC AEZ2 A CH20 – DRC AEZ2 B CH23 Uganda (H8: KX276941, KU725897) – B CH24 Kenya (H14: KU725903) – B CH25 Kenya (H27: KU725916) – B f,j e CH26 Grande Comore (H2: DQ901357) , Kenya (H7: KU725896, H13: KU725902) – B CH27 Kenya (H16: KU725905) – B e f f CH28 Kenya (H21: KU725910 , H9: DQ901359 , H9: DQ901358 ) – B CH29 Kenya (H28: KU725917) – B CH30 Kenya (H9: KU725898) – A CH31 Kenya (H17: KU725906 – A CH32 Kenya (H24: KU725913) – A CH33 – Tanzania (TZ13) A CH34 Kenya (H22: KU725911) – A CH35 Kenya (H23: KU725912) – A CH36 Kenya (H10: KU725899) – A CH37 Kenya (H15:KU725904) – A CH38 Kenya (H20: KU725909) – A CH39 Kenya (H19: KU725908) – A CH40 Kenya (H26: KU725915) – A CH41 Kenya (H25: KU725914 – A CH42 Kenya (H18: KU725907) – A CH43 Kenya (H12: KU725901) – A CH44 Zambia-east (H3: DQ859263) – A CH45 Zambia-east (H4: DQ859264) – A CH1: variants CH1, CH2, CH5, CH15, CH17 and CH21 (Additional file 3: Table S3) CH6: variants CH6, CH8 and CH22 (Additional file 3: Table S3) CH7: variants CH7 and CH13; d CH11: CH11 and CH19 (Additional file 3: Table S3) Haplotypes underlined are exclusive to the Great Lakes region. Similar data for 12S rRNA are detailed in Additional file 7: Table S6 e f g h i j [31]; [29]; [52]; [68]; [69]; [14] Amzati et al. Parasites & Vectors (2018) 11:329 Page 15 of 19 continuous rainfall rather harbour smaller ticks without tick R. appendiculatus that took place in the Great Lakes diapause with bivoltine or continuous phenology, corre- region. A strong evidence of recent spatial and demo- sponding to “east African group” or haplogroup A [37, 52]. graphic expansion was observed for lineage A in all AEZs The highest genetic variability was observed in low- included in the study. Analyses of cox1 sequences revealed lands, whereas a relatively lower diversity was observed relatively high haplotype diversity contrasted with low nu- in midlands and highlands. The high genetic diversity in cleotide diversity values for each population, suggesting a lowlands can be explained by the strong presence of the sudden population expansion [54, 55]. This result, together two lineages A and B observed in these AEZs. The coex- with negative values of Tajima’s D and Fu’s Fs, the star-like istence of these lineages could originate from the disper- radiation, the unimodal mismatch pattern and sal of the tick through livestock movement between non-significant RI and SSD statistics, further support the re- AEZs [53], associated to the suitability of semi-arid cli- cent evolutionary history andsuddenpopulationgrowth mate for lineage B expressing obligatory diapause [37]. experienced by lineage A [56–58]. We did not estimate the expansion time, but we hypothesize that the expansion was Moderate genetic structure of R. appendiculatus between recent and not sufficient to increase the nucleotide diver- lowlands and highlands sity, probably because of recent coalescence, while rapid Population genetic analyses of cox1 gene variation in R. population expansion following a selective sweep (bottle- appendiculatus revealed low to moderate genetic differen- neck or genetic drift) could have accumulated new muta- tiation values and high gene flow rates among AEZs. The tions that sufficiently increased the haplotype diversity [45, two identified R. appendiculatus lineages were sympatric 59, 60]. This could explain the excess of singletons poly- in the Great Lakes region, although lineage A was the morphic sites and rare haplotypes that diverge from most abundant and widely distributed in all AEZs and ancestral haplotypes by only 1–2mutationalsteps [45, 61]. lineage B was particular confined to lowlands, where the However, the strong bimodal mismatch pattern observed in climate is tropical dry, more arid with lower annual rain- Rwanda AEZ2 and DRC AEZ1 suggests the coexistence of fall and longer dry season than in highlands (short dry sea- R. appendiculatus lineages after recent colonization events son with abundant annual rainfall). These climate or exchanging migrants [45, 56, 62]. The scenario of conditions in lowlands are quite similar to the described lineage B contrasts with that of lineage A. Analyses indi- ecological zones of lineage B in southern Africa [29, 52]. cate, for lineage B, an equilibrium situation that is charac- In addition, altitudinal gradient seems to be a key factor terised by a weak signal of recent bottleneck and no that shapes the distribution pattern and the establishment evidence of population expansion. It was also less diverse ability of lineage B. Its frequency decreases with increasing than haplogroup A, indicating that haplogroup A is ex- altitude. High degree of genetic similarity was observed periencing population expansion independently of hap- between the lowlands of DRC and the low plateau of logroup B and it has been established longer in the Great Rwanda and between the highlands of DRC and Burundi, Lakes region, while haplogroup B was probably which are geographically distant from each other. The introduced more recently and established a founder popu- most likely explanation for this is that the spatial pattern lation. There are three possible explanations of the equi- of R. appendiculatus lineages is not only driven by geo- librium observed in haplogroup B: (i) only few haplotypes graphical separation as described in previous studies [52], were recently introduced; (ii) only the few identified hap- but also related to their ecological preferences, as ob- lotypes persisted after a bottleneck; or (iii) haplogroup B is served by the significant genetic differentiation among experiencing an initial selection sweep which has reduced lowlands and highlands AEZs. On the other hand, adja- the number of rare haplotypes and singleton mutations cent AEZs shared more migrants, especially of lineage A, [56]. When we analysed R. appendiculatus cox1sequences facilitated by short-distance seasonal movement of cattle available in Africa, the star radiation of the MJ network in [9], which may have reduced the geographical structuring each group suggests that the two lineages went through a of the tick [11, 18]. Analysis of molecular variance demographic expansion and evolve independently of each (AMOVA) confirmed these findings showing that the vari- other with limited gene exchange. Unfortunately, we were ance explained by divergence between the six AEZs was not able to test the hypothesis of crossbreeding events be- lower (6%), while the largest fraction of genetic variation tween the two lineages described here, because of the ma- was observed among individuals within AEZs (94%). ternal inheritance of mitochondrial genes used in the present study. More studies, such as extensive biological Rhipicephalus appendiculatus lineage A has undergone characterisation, crossbreeding experiments and the use sudden demographic and range expansion in the Great of biparental inheritance markers are necessary to investi- Lakes region gate the panmixia of the two lineages and the genetic The demographic and spatial dynamics were analysed using mechanisms driving their establishment and correspond- multiple algorithms, to elucidate colonization events of the ing phenotypic variations in changed environments. Amzati et al. Parasites & Vectors (2018) 11:329 Page 16 of 19 The fact that R. appendiculatus has undergone spatial behaviour for southern African ticks, which allow them expansion was in accordance with the expectation that to survive hot and dry ecological conditions [33, 37]. We long and short distance movement of cattle are key factors hypothesise that these characteristics make the develop- of spreading ticks. The processes responsible for that evo- ment of lineage B slower than lineage A in sympatric lutionary pattern may have resulted, not only from range areas, giving an evolutionary advantage to the latter. expanding of the previously established haplotypes to This could also reduce the abundance of lineage B delay- proximate AEZs, but also from the recolonization events ing its oviposition during unfavourable conditions [37]. of ticks from other regions and countries [53, 63]. In The processes that took place to divide the two groups addition to cattle movement, the population expansion need to be further investigated. However, we propose and establishment of novel haplotypes in previously un- that they could have diverged following genetic drift due occupied areas could have been driven by recent environ- to founder events of natural geographical barrier that mental and climatic changes, affecting vector-borne may result in reproductive isolation. It is demonstrated diseases landscape over recent decades [64–66]. that biogeographical break again host migration reduces gene exchange and could dictate the spatial and repro- Sympatric and allopatric ecological zones of R. ductive separation within a species [60]. appendiculatus lineages in Africa Phylogenetic trees produced two main genetic subpopu- Conclusions lations of R. appendiculatus that have a wide distribution This study provided new insights into the genetic struc- range in Africa, with large divergence in behavioural dia- ture and colonization events of R. appendiculatus in the pause [36, 37], spatial variation in body size [35, 51], Great Lakes region and over its distribution range in ecology and phenology [32–34]. Initially, the two line- sub-Saharan Africa. Our results highlighted the occur- ages were considered as “east African” and “southern Af- rence of two major genetic lineages (A and B) in the Great rican” genetic groups [29]. To date, they are sympatric Lakes region. The two lineages are not spatially structured in most eastern and central African ecological zones. For in the study region and they differ in their colonization instance, previous studies did not reveal the presence of histories and pattern. Lineage B, not previously reported, haplogroup B in Rwanda [29]. This could be an indica- was probably introduced recently in the region and its oc- tion of recent introduction of the tick in the Great Lakes currence decreases with increased altitude, whereas region. The MJ network further elucidate that the initial lineage A, widely distributed, has been longer established population of haplogroup B could have come from an and subjected to sudden demographic and spatial expan- ancestral female of haplotype CH7, which is the most sion most likely related to short and long-distance cattle prevalent haplotype of this group in areas where it oc- movement. Rhipicephalus appendiculatus ticks are more curs. Consequently, two different eco-genetic zones are diverse in lowlands than highlands with moderate genetic shaped in Africa, the sympatric zone where the two line- differentiation between the two ecosystems, while more ages are found, which covers central and eastern Africa, genetic similitude is found in zones with same and allopatric zone in southern Africa where the two lin- agro-ecological profiles, in spite of their geographical dis- eages have clear geographical and ecological separation. tance. The genetic distribution of R. appendiculatus sug- For instance, in Zambia, lineage B is present in the south gests two different eco-genetic zones in Africa, the (long dry season) and lineage A in the east of the coun- sympatric zone (central and eastern Africa) where the two try (shorter dry season) [29]. In Grande Comore, lineage lineages coexist and the allopatric zones (southern Africa) B has established a stable population, while lineage A where they have clear geographical divergence. The range was found on imported cattle [14]. In areas where the expansion pattern of lineages and the genetic admixture two lineages are sympatric, their respective abundances of R. appendiculatus populations observed in the Great differ, mainly driven by their divergence in ecological Lakes region can strongly affect the epidemiological dy- preferences and plasticity [33]. The sympatric relation- namics of ECF. This could partially explain the endemic ship is in agreement with the observation made by Berk- instability and occasional epidemics due to the introduc- vens et al. [67], where an east African stock from Kenya tion and temporal subsistence of infected ticks mostly in expressed diapause, contrasting with the result obtained fringes areas of lowlands. by Madder et al. [36], where another stock from the same region was unable to enter diapause. These evi- Additional files dences show that lineage B has a greater invasive ability into new habitats and better fit wide range of tropical Additional file 1: Table S1. Rhipicephalus appendiculatus cox1 and 12S rRNA haplotype sequences retrieved from GenBank. (DOCX 16 kb) and equatorial conditions, while lineage A is particularly Additional file 2: Table S2. cox1 and 12S rRNA BLAST results for species confined to equatorial conditions. This could be ex- identification and confirmation. (DOCX 16 kb) plained by the larger body size and obligatory diapause Amzati et al. Parasites & Vectors (2018) 11:329 Page 17 of 19 GenBank under accession numbers MF458950-MF458971 and the GenBank Additional file 3: Table S3. Polymorphism in the 22 haplotypes of the numbers for the nine 12S rRNA haplotypes are MF479189-MF479197. cox1 gene fragment of R. appendiculatus. (DOCX 21 kb) Additional file 4: Table S4. Population genetic structure inferred by Authors’ contributions analysis of molecular variance (AMOVA) based on cox1 sequences of R. GSA collected and identified ticks, performed experiments, analysed data appendiculatus from different agro-ecological zones. (DOCX 13 kb) and wrote the manuscript. TM, NK and GSA conceived the study. TM, NK, RP, JBM, AD and MM participated in the study design, supervised the research Additional file 5: Table S5. Evolutionary neutrality, demographic and and critically revised the manuscript. RP and EGK participated in laboratory spatial history of mitochondrial cox1 gene. (DOCX 15 kb) experiments and advised in data analyses. JBM supervised tick collections. All Additional file 6: Figure S1. cox1 mismatch distribution pattern for R. authors read and approved the final manuscript. appendiculatus haplogroup A in different agro-ecological zones. (DOCX 193 kb) Additional file 7: Table S6. Rhipicephalus appendiculatus 12S rRNA Ethics approval and consent to participate haplotypes and their distribution among agro-ecological zones of the The research was carried out in accordance with ethical guidelines and Great lakes region and other sub-Saharan African countries. (DOCX 15 kb) regulations of the Université Evangélique en Afrique (UEA/Bukavu-DR Congo). The specific field protocols were approved by the UEA (Reference Additional file 8: Figure S2. Neighbor-joining tree of 12S haplotype number: UEA/SGAC/KM/020/2015) and permitted by the Provincial sequences for R. appendiculatus across African countries. (DOCX 18 kb) Inspection of Agriculture, Livestock and Fisheries (Reference number: 55.00/ 0026/IPAPEL/SK/2015). Farm managers and owners were informed about the Abbreviations survey and gave approval for tick collection on their cattle. This research did AMOVA: Analysis of molecular variance; BecA: Biosciences eastern and central not involve species at risk of extinction. Africa; BLAST: Basic local alignment search tool; cox1: Cytochrome c oxidase subunit 1; D: Tajima’s neutrality test; ECF: East Coast fever; Fs:Fu’s neutrality Competing interests tests; F : Right’s fixation index; h: Number of haplotypes; Hd: Haplotype ST The authors declare that they have no competing interests. diversity; ILRI: International Livestock Research Institute; ITS2: Internal transcribed spacer 2; K: Mean number of pairwise nucleotide differences Publisher’sNote within population; K : Average number of nucleotide differences between XY Springer Nature remains neutral with regard to jurisdictional claims in published populations; ML: Maximum likelihood; NCBI: National Center for maps and institutional affiliations. Biotechnology Information; NJ: Neighbor-joining; Nm: Number of migrants; PCR: Polymerase chain reaction; PIS: Parsimony informative sites; Author details RI: Harpending’s raggedness index; rRNA: ribosomal RNA; S: Segregating sites; Unit of Integrated Veterinary Research, Department of Veterinary Medicine, SD: Standard deviation; SSD: Sum of squares deviation; π: nucleotide diversity Faculty of Sciences, Namur Research Institute for Life Sciences (NARILIS), University of Namur (UNamur), Rue de Bruxelles 61, 5000 Namur, Belgium. Acknowledgments Research Unit of Veterinary Epidemiology and Biostatistics, Faculty of We would like to thank Professor Paul Gwakisa from the Sokoine University Agricultural and Environmental Sciences, Université Evangélique en Afrique, of Agriculture (SUA) in Tanzania for providing tick specimens from Tanzania. P.O. Box 3323, Bukavu, Democratic Republic of the Congo. Biosciences We thank Alphone Bisusa from the Centre de Recherche en Sciences eastern and central Africa - International Livestock Research Institute Naturelles (CRSN/Lwiro) and Amzati Saidi Ngton (Université Evangélique en (BecA-ILRI) hub, P.O. Box 30709-00100, Nairobi, Kenya. Department of Afrique) for technical assistance with sample collection and tick Biochemistry, School of Medicine, University of Nairobi, P.O. Box identification. We are also grateful to the technical staff of BecA-ILRI for facili- 30197-00100, Nairobi, Kenya. Present address: Centre for Tropical Livestock tating laboratory work, especially Stephen Mwaura (tick unit-ILRI) for assist- Genetics and Health (CTLGH), The University of Edinburgh, Easter Bush, ance with the confirmation of morphological identification of closely related Midlothian, Scotland EH25 9RG, UK. Department of Veterinary Tropical tick species. We thank Dr Marie Cariou (Research Unit in Environmental and Diseases, Faculty of Veterinary Science, University of Pretoria, P/Bag X04, Evolutionary Biology, University of Namur) for advising on data analysis and Onderstepoort 0110, South Africa. interpretation. We are also grateful to Professor Gustave Mushagalusa (Uni- versité Evangélique en Afrique) for mobilizing partnerships and collabora- Received: 11 January 2018 Accepted: 16 May 2018 tions to the “Theileria” project in DRC. Funding References The laboratory aspects of this project were fully supported by the BecA-ILRI 1. Nene V, Kiara H, Lacasta A, Pelle R, Svitek N, Steinaa L. The biology of Hub through the Africa Biosciences Challenge Fund (ABCF) program. The Theileria parva and control of East Coast fever - Current status and future ABCF Program is funded by the Australian Department for Foreign Affairs trends. Ticks Tick Borne Dis. 2016;7:549–64. and Trade (DFAT) through the BecA-CSIRO partnership; the Syngenta Foun- 2. Kalume MK, Saegerman C, Mbahikyavolo DK, Makumyaviri AM, Marcotty T, dation for Sustainable Agriculture (SFSA); the Bill & Melinda Gates Foundation Madder M, et al. Identification of hard ticks (Acari: Ixodidae) and (BMGF); the UK Department for International Development (DFID) and; the seroprevalence to Theileria parva in cattle raised in North Kivu Province, Swedish International Development Cooperation Agency (Sida). Sample col- Democratic Republic of Congo. Parasitol Res. 2013;112:789–97. lection, field equipment and scientific workshop organised in DRC were sup- 3. Bazarusanga T, Geysen D, Vercruysse J, Madder M. An update on the ported through the “Theileria” project co-funded to the Université ecological distribution of ixodid ticks infesting cattle in Rwanda: Evangélique en Afrique (UEA) by the Agence Universitaire de la Francopho- countrywide cross-sectional survey in the wet and the dry season. Exp Appl nie (AUF) and the Communauté Economique des Pays des Grands Lacs Acarol. 2007;43:279–91. (CEPGL). The International foundation of Science (IFS) supported in collabor- 4. Kaiser M, Sutherst R, Bourne A, Gorissen L, Floyd R. Population dynamics of ation with BecA-ILRI hub the individual scholarship awarded to GSA for field ticks on Ankole cattle in five ecological zones in Burundi and strategies for work and part of field equipment to the “Theileria” project. The study was their control. Prev Vet Med. 1988;6:199–222. also supported by the University of Namur through the UNamur-CERUNA in- 5. Olwoch J, Reyers B, van Jaarsveld A. Host-parasite distribution patterns stitutional PhD grant awarded to GSA for manuscript write-up in Belgium. under simulated climate: implications for tick-borne diseases. Int J Climatol. 2009;29:993–1000. Availability of data and materials 6. Perry B, Lessard P, Norval R, Kundert K, Kruska R. Climate, vegetation The data supporting the findings of this study are included within the article and the distribution of Rhipicephalus appendiculatus in Africa. Parasitol and its additional files. Representative sequences of R. appendiculatus from Today. 1990;6:100–4. each agro-ecological zone (AEZ) are available in the GenBank database under 7. Mararo SB. Pouvoirs, élevage bovin et la question foncière au Nord-Kivu. In: accession numbers: MF458895-MF458949 and MF479166-MF479188 for cox1 Marysse S, Reyntjens F, editors. L'Afrique des Grands Lacs: annuaire 2000– and 12S rRNA genes, respectively. The 22 cox1 haplotypes were submitted to 2001. Paris: L'Harmattan; 2001. p. 219–50. Amzati et al. Parasites & Vectors (2018) 11:329 Page 18 of 19 8. De Failly D. L’économie du Sud-Kivu 1990–2000: mutations profondes 33. Speybroeck N, Madder M, Van Den Bossche P, Mtambo J, Berkvens N, Chaka cachées par une panne. In: Marysse S, Reyntjens F, editors. L’Afrique des G, et al. Distribution and phenology of ixodid ticks in southern Zambia. Med Grands Lacs: annuaire 1999–2000. Paris: L'Harmattan; 2000. p. 163–92. Vet Entomol. 2002;16:430–41. 9. Verweijen J, Brabant J. Cows and guns. Cattle-related conflict and 34. Leta S, De Clercq EM, Madder M. High-resolution predictive mapping for armed violence in Fizi and Itombwe, eastern DR Congo. J Mod Afr Rhipicephalus appendiculatus (Acari: Ixodidae) in the Horn of Africa. Exp Appl Stud. 2017;55:1–27. Acarol. 2013;60:531–42. 10. Vlassenroot K, Huggins C. Land, migration and conflict ineasternDRC. In: 35. Chaka G, Billiouw M, Geysen DM, Berkvens DL. Spatial and temporal Huggins C, Clover J, editors. From the ground up: land rights, conflict and peace variation in Rhipicephalus appendiculatus size in eastern Zambia. Trop Med in sub-Saharan Africa. Pretoria: Institute for Security Studies; 2005. p. 115–94. Int Health. 1999;4:A43–8. 11. Boulinier T, Kada S, Ponchon A, Dupraz M, Dietrich M, Gamble A, et al. 36. Madder M, Speybroeck N, Brandt J, Berkvens D. Diapause induction in Migration, prospecting, dispersal? What host movement matters for adults of three Rhipicephalus appendiculatus stocks. Exp Appl Acarol. 1999; infectious agent circulation? Integr Comp Biol. 2016;56:330–42. 23:961–8. 12. Maze-Guilmo E, Blanchet S, McCoy KD, Loot G. Host dispersal as the driver 37. Madder M, Speybroeck N, Brandt J, Tirry L, Hodek I, Berkvens D. Geographic of parasite genetic structure: a paradigm lost? Ecol Lett. 2016;19:336–47. variation in diapause response of adult Rhipicephalus appendiculatus ticks. 13. De Deken R, Martin V, Saido A, Madder M, Brandt J, Geysen D. An outbreak Exp Appl Acarol. 2002;27:209–21. of East Coast fever on the Comoros: a consequence of the import of 38. Bazarusanga T, Vercruysse J, Marcotty T, Geysen D. Epidemiological studies immunised cattle from Tanzania? Vet Parasitol. 2007;143:245–53. on theileriosis and the dynamics of Theileria parva infections in Rwanda. Vet 14. Yssouf A, Lagadec E, Bakari A, Foray C, Stachurski F, Cardinale E, et al. Parasitol. 2007;143:214–21. Colonization of Grande Comore Island by a lineage of Rhipicephalus 39. Walker AR, Bouattour A, Camicas J-L, Estrada-Pena A, Horak IG, Latif AA, appendiculatus ticks. Parasit Vectors. 2011;4:38. et al. Ticks of domestic animals in Africa: a guide to identification of species. Edinburgh, Scotland: Bioscience Reports; 2003. 15. Fevre EM, Bronsvoort BM, Hamilton KA, Cleaveland S. Animal movements and the spread of infectious diseases. Trends Microbiol. 2006;14:125–31. 40. Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R. DNA primers for 16. Barre N, Uilenberg G. Spread of parasites transported with their hosts: case amplification of mitochondrial cytochrome c oxidase subunit I from diverse study of two species of cattle tick. Rev Sci Tech. 2010;29:135–47. metazoan invertebrates. Mol Mar Biol Biotechnol. 1994;3:294–9. 17. Estrada-Peña A, Salman M. Current limitations in the control and spread of 41. Simon C, Frati F, Beckenbach A, Crespi B, Liu H, Flook P. Evolution, ticks that affect livestock: a review. Agriculture. 2013;3:221–35. weighting, and phylogenetic utility of mitochondrial gene sequences and a 18. Leo SS, Gonzalez A, Millien V. The genetic signature of range expansion in a compilation of conserved polymerase chain reaction primers. Ann Entomol disease vector-the black-legged tick. J Hered. 2017;108:176–83. Soc Am. 1994;87:651–701. 42. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA 19. Ogden N. Changing geographic ranges of ticks and tick-borne pathogens: polymorphism data. Bioinformatics. 2009;25:1451–2. drivers, mechanisms and consequences for pathogen diversity. Front Cell 43. Excoffier L, Lischer HE. Arlequin suite ver 3.5: a new series of programs to Infect Microbiol. 2013;3:46. perform population genetics analyses under Linux and Windows. Mol Ecol 20. Criscione CD, Poulin R, Blouin MS. Molecular ecology of parasites: elucidating Resour. 2010;10:564–7. ecological and microevolutionary processes. Mol Ecol. 2005;14:2247–57. 21. Gooding RH. Genetic variation in arthropod vectors of disease-causing 44. Meirmans PG, Hedrick PW. Assessing population structure: FST and related organisms: obstacles and opportunities. Clin Microbiol Rev. 1996;9:301–20. measures. Mol Ecol Resour. 2011;11:5–18. 22. Gandon S, Michalakis Y. Local adaptation, evolutionary potential and host- 45. Rogers AR, Harpending H. Population growth makes waves in the parasite coevolution: interactions between migration, mutation, population distribution of pairwise genetic differences. Mol Biol Evol. 1992;9:552–69. size and generation time. J Evol Biol. 2002;15:451–62. 46. Harpending HC. Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Hum Biol. 1994;66:591–600. 23. Kubasu S. The ability of Rhipicephalus appendiculatus (Acarina: Ixodidae) 47. Fu YX. Statistical tests of neutrality of mutations against population growth, stocks in Kenya to become infected with Theileria parva. Bull Entomol Res. hitchhiking and background selection. Genetics. 1997;147:915–25. 1992;82:349–53. 48. Tajima F. Statistical method for testing the neutral mutation hypothesis by 24. Odongo DO, Ueti MW, Mwaura SN, Knowles DP, Bishop RP, Scoles GA. DNA polymorphism. Genetics. 1989;123:585–95. Quantification of Theileria parva in Rhipicephalus appendiculatus (Acari: Ixodidae) confirms differences in infection between selected tick strains. J 49. Kumar S, Stecher G, Tamura K. MEGA7: Molecular Evolutionary Genetics Med Entomol. 2009;46:888–94. Analysis Version 7.0 for bigger datasets. Mol Biol Evol. 2016;33:1870–4. 25. Young AS, Dolan TT, Mwakima FN, Ochanda H, Mwaura SN, Njihia GM, et al. 50. Leigh JW, Bryant D. popart: full-feature software for haplotype network Estimation of heritability of susceptibility to infection with Theileria parva in construction. Methods Ecol Evol. 2015;6:1110–6. the tick Rhipicephalus appendiculatus. Parasitology. 1995;111:31–8. 51. Speybroeck N, Madder M, Thulke HH, Mtambo J, Tirry L, Chaka G, et al. 26. Ochanda H, Young A, Medley G, Perry B. Vector competence of 7 Variation in body size in the tick complex Rhipicephalus appendiculatus/ rhipicephalid tick stocks in transmitting 2 Theileria parva parasite stocks Rhipicephalus zambeziensis. J Vector Ecol. 2004;29:347–54. from Kenya and Zimbabwe. Parasitology. 1998;116:539–45. 52. Mtambo J, Madder M, Van Bortel W, Chaka G, Berkvens D, Backeljau T. Further evidence for geographic differentiation in R. appendiculatus (Acari: 27. McCoy K. The population genetic structure of vectors and our Ixodidae) from Eastern and Southern provinces of Zambia. Exp Appl Acarol. understanding of disease epidemiology. Parasite. 2008;15:444–8. 2007;41:129–38. 28. Le Roux J, Wieczorek A. Molecular systematics and population genetics of biological invasions: towards a better understanding of invasive species 53. Excoffier L, Foll M, Petit RJ. Genetic consequences of range expansions. management. Ann Appl Biol. 2009;154:1–7. Annu Rev Ecol Evol Syst. 2009;40:481–501. 29. Mtambo J, Madder M, Van Bortel W, Geysen D, Berkvens D, Backeljau T. 54. Braverman JM, Hudson RR, Kaplan NL, Langley CH, Stephan W. The Genetic variation in Rhipicephalus appendiculatus (Acari: Ixodidae) from Zambia: hitchhiking effect on the site frequency spectrum of DNA polymorphisms. correlating genetic and ecological variation with Rhipicephalus appendiculatus Genetics. 1995;140:783–96. from eastern and southern Africa. J Vector Ecol. 2007;32:168–75. 55. Simonsen KL, Churchill GA, Aquadro CF. Properties of statistical tests of 30. Kanduma EG, Mwacharo JM, Mwaura S, Njuguna JN, Nzuki I, Kinyanjui PW, neutrality for DNA polymorphism data. Genetics. 1995;141:413–29. et al. Multi-locus genotyping reveals absence of genetic structure in field 56. Ray N, Currat M, Excoffier L. Intra-deme molecular diversity in spatially populations of the brown ear tick (Rhipicephalus appendiculatus) in Kenya. expanding populations. Mol Biol Evol. 2003;20:76–86. Ticks Tick Borne Dis. 2016;7:26–35. 57. Frankham R. Relationship of genetic variation to population size in wildlife. Conserv Biol. 1996;10:1500–8. 31. Kanduma EG, Mwacharo JM, Githaka NW, Kinyanjui PW, Njuguna JN, Kamau LM, et al. Analyses of mitochondrial genes reveal two sympatric but 58. Emerson BC, Paradis E, Thébaud C. Revealing the demographic histories of genetically divergent lineages of Rhipicephalus appendiculatus in Kenya. species using DNA sequences. Trends Ecol Evol. 2001;16:707–16. Parasit Vectors. 2016;9:353. 59. Robbertse L, Baron S, van der Merwe NA, Madder M, Stoltsz WH, Maritz- 32. Berkvens DL, Geysen DM, Chaka G, Madder M, Brandt JR. A survey of the Olivier C. Genetic diversity, acaricide resistance status and evolutionary ixodid ticks parasitising cattle in the Eastern Province of Zambia. Med Vet potential of a Rhipicephalus microplus population from a disease-controlled Entomol. 1998;12:234–40. cattle farming area in South Africa. Ticks Tick Borne Dis. 2016;7:595–603. Amzati et al. Parasites & Vectors (2018) 11:329 Page 19 of 19 60. Avise JC. Phylogeography: The history and formation of species. USA: Harvard University Press; 2000. 61. Slatkin M, Hudson RR. Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations. Genetics. 1991; 129:555–62. 62. Rogers AR, Fraley AE, Bamshad MJ, Watkins WS, Jorde LB. Mitochondrial mismatch analysis is insensitive to the mutational process. Mol Biol Evol. 1996;13:895–902. 63. Cristescu ME. Genetic reconstructions of invasion history. Mol Ecol. 2015;24:2212–25. 64. Biek R, Real LA. The landscape genetics of infectious disease emergence and spread. Mol Ecol. 2010;19:3515–31. 65. Ostfeld RS. Climate change and the distribution and intensity of infectious diseases. Ecology. 2009;90:903–5. 66. Manel S, Holderegger R. Ten years of landscape genetics. Trends Ecol Evol. 2013;28:614–21. 67. Berkvens DL, Pegram RG, Brandt JR. A study of the diapausing behaviour of Rhipicephalus appendiculatus and R. zambeziensis under quasi-natural conditions in Zambia. Med Vet Entomol. 1995;9:307–15. 68. Murrell A, Campbell NJH, Barker SC. Phylogenetic analyses of the rhipicephaline ticks indicate that the genus Rhipicephalus is paraphyletic. Mol Phylogenet Evol. 2000;16:1–7. 69. Burger TD, Shao R, Barker SC. Phylogenetic analysis of mitochondrial genome sequences indicates that the cattle tick, Rhipicephalus (Boophilus) microplus, contains a cryptic species. Mol Phylogenet Evol. 2014;76:241–53.

Journal

Parasites & VectorsSpringer Journals

Published: May 31, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off