Accurate Phylogenetic Relationships Among Mycobacterium bovis Strains Circulating in France Based on Whole Genome Sequencing and Single Nucleotide Polymorphism Analysis

In recent years the diversity of the French Mycobacterium bovis population responsible for bovine tuberculosis (bTB) outbreaks since 1970 has been described in detail. To further understand bTB evolution in France, we used single nucleotide polymorphisms (SNPs) based on whole genome sequence versus classical genotyping methods in order to identify accurate phylogenetic relationships between M. bovis strains. Whole genome sequencing was carried out on a selection of 87 strains which reflect the French M. bovis population’s genetic diversity. Sequences were compared to the M. bovis reference genome AF2122/97. Comparison among the 87 genomes revealed 9,170 sites where at least one strain shows a SNP with respect to the reference genome; 1,172 are intergenic and 7,998 in coding sequences, of which 2,880 are synonymous and 5,118 non-synonymous. SNP-based phylogenetic analysis using these 9,170 SNP is congruent with the cluster defined by spoligotyping and multilocus variable number of tandem repeat analysis typing. In addition, some SNPs were identified as specific to genotypic groups. These findings suggest new SNP targets that can be used for the development of high-resolving methods for genotyping as well as for studying M. bovis evolution and transmission patterns. The detection of non-synonymous SNPs on virulence genes enabled us to distinguish different clusters. Our results seem to indicate that genetically differentiated clusters could also display distinctive phenotypic traits.


INTRODUCTION
Mycobacterium bovis, the main etiological agent of bovine tuberculosis (bTB), is a member of the Mycobacterium tuberculosis complex (MTBC) (Huard et al., 2006;Smith et al., 2006a,b). This highly clonal species, which derived from a common ancestor related to Mycobacterium canettii (Blouin et al., 2012), includes seven different members: M. tuberculosis, M. bovis, M. africanum, M. microti, M. caprae, M. canettii, M. pinnipedii. Members of the M. tuberculosis complex may be considered a series of host-adapted ecotypes rather than species that are part of a strain lineage (Smith et al., 2006a(Smith et al., , 2009. These ecotypes are characterized by the absence of a chromosomal region. M. bovis is the final ecotype of a series of evolution from an ancestor related to M. tuberculosis. To date, four clonal complexes of M. bovis have been defined on the basis of deletions, distinct spoligotype and SNPs signatures: European 1, European 2, African 1, African 2 (Muller et al., 2009;Smith, 2012;Rodriguez-Campos et al., 2014). In France, as in many countries, bTB has a great economic impact due to the high cost of the control and eradication campaign as well as the impossibility to trade live animals or animal products (Anonymous, 2000;Amanfu, 2006). Furthermore, M. bovis infection is a zoonosis and a major concern at the animal-human interface, mainly through the consumption of raw milk or raw milk products (Palmer et al., 2012;Muller et al., 2013) particularly in those regions of the world where predisposing factors such as poverty or aids prevail.
Genetic differentiation of strains has become an indispensable tool to study the evolution, epidemiology and ecology of this pathogen.
For epidemiological surveillance, spoligotyping and the multilocus variable number of tandem repeat analysis (MLVA), two molecular tools that make it possible to discriminate strains by their genotypes (Aranaz et al., 1996;Skuce et al., 2005;Allix et al., 2006;Supply et al., 2006) are currently used worldwide (Haddad et al., 2001;Boniotti et al., 2009;Rodriguez et al., 2010;Smith, 2012;Rodriguez-Campos et al., 2013). Using these conventional methods, we recently described the dynamic genetic evolution of the French M. bovis population studying strains isolated mainly from livestock and wildlife that has been isolated since 1978 (Hauer et al., 2015). During this period, the number of bTB-positive herds decreased from 15% to less than 0.1% although outbreaks persisted in some regions. This decrease in prevalence was concomitant with a decrease of the genetic diversity of strains characterized by a strong regionalization and the evidence of multi-host transmission cycles in which both livestock and wildlife are infected by strains of the same genotype. The observed genetic relationships between the strains supports the co-existence of at least four clonal populations within the French M. bovis population (Hauer et al., 2015).
Whilst spoligotyping and MLVA are unquestioned tools for disease outbreak investigations, recent studies show the limitations of these methods for tracing the origin of a given infection (Rodriguez-Campos et al., 2011). Although these typing techniques target polymorphic genetic regions, they interrogate less than 1% of the genome and therefore have an intrinsically restricted discriminatory power. Furthermore regarding spoligotyping, the deletion of the same type of spacer can occur in entirely separate sub-lineages (Corner et al., 1995). Likewise, copy numbers of MIRU (Mycobacterial Interspersed Repetitive Units) -VNTR (Variable Numbers of Tandem Repeats) at each locus can increase or decrease during DNA replication which may contribute to give rise to identical patterns in separate sub-lineages (Corner et al., 1995).
These limitations can be overcome by the application of whole genome sequencing (WGS) for genome-based epidemiology especially for clonal populations belonging to highly concordant phylogenies with low homoplasy. WGS can provide comprehensive genetic information including all possible genomic targets as well as additional information on genome evolution, virulent determinants and drug resistance traits. Using the WGS, single nucleotide polymorphisms (SNPs) analyses become a robust tool for studying closely related strains of pathogenic mycobacteria Homolka et al., 2012;Joshi et al., 2012;Roetzer et al., 2013).
The purpose of the present study was to undertake SNP analysis based on WGS versus classic genotyping for identifying accurate phylogenetic relationships within the population of French M. bovis strains. A set of 87 animal field strains representative of the French M. bovis population's genetic diversity was selected based on a previous large-scale genotyping study (Hauer et al., 2015). We performed WGS of the 87 field isolates to detect SNPs by genome comparison among these M. bovis strains using the genome M. bovis AF2122/97 as reference. The SNP analysis was addressed to verify the accuracy of the phylogenetic lineages of the isolates and to establish correlations between genotypes and phenotypes based on particular SNP profiles targeting virulence genes.  (Table 1 and  Supplementary Table S1) was carefully selected so as to represent the genetic diversity of the main genotypes described in our recent previous study. Isolates were obtained from animals studied in the framework of French bTB control campaign from different hosts: wild boars (n = 2); pig (n = 1); sheep (n = 1); goat (n = 1); deer (n = 1) and mostly from cattle (n = 81). Details of the spoligotype and MLVA profiles of the isolates selected to be sequenced are indicated in Table 1 and Supplementary Figure S1. The following six query genomes were also included as references: M. bovis AF2122/97_BX24833 (EMBL sequence AF2122/97_BX248333); M. tuberculosis H37Rv (NC_000962.2); M. africanum MAL010084 (JKXR00000000.1); M. canettii CIPT 140010059 (NC_015848.1); M. bovis BCG Pasteur_1173P2 (NC_008769.1); and M. bovis BCG Tokyo_172 (NC_012207.1).

DNA Extraction
Clones of the different isolates were enriched in 50 mL Middlebrook 7H9 broth and harvested at mid log phase growth and pelleted at room temperature at 5,500 g. Pellets were inactivated by heating at 80 • C during 30 min. DNA was extracted by a chloroform protocol with a lysozyme pre-step followed by a proteinase K treatment. CTAB/NACL was added before the chloroform/isoamyl alcohol action. Precipitation of DNA was performed in isopropanol. DNA was conserved at −20 • C in 1x TE buffer (Imai et al., 1994;Biek et al., 2012).

Genotyping
Spoligotyping of 43 spacers sequences contained in the direct repeat (DR) locus was performed on a Luminex (Zhang et al., 2010). Spoligotypes were named according to an agreed international convention 1 (Smith, 2012). MLVA typing (Frothingham and Meeker-O'Connell, 1998;Roring et al., 2002) was performed by Genoscreen (Lille, France), on 8 MIRU -  Table S2). To avoid PCR overrepresented fragments during the library preparation, the paired-end FASTQ files were filtered from replicated sequences, leaving only one pair of short reads for each set of pair of replicated short reads. The short reads were trimmed using Trimmomatic (Bolger et al., 2014) when the quality call given by the sequencer (range: 0 to 40) across the sequence dropped lower than 20 for a slidding window of 10bp (Bolger et al., 2014). The minimum read length threshold equal 35 bases. The resulting mean length of the short reads per sample ranged from 95.31 to 98.42 bp with standard deviation ranging from 10.25 to 13.78 (Supplementary Table S2).
The FASTQ were then mapped against the M. bovis reference strain AF2122/97-BX24833 using SMALT (Sanger Institute). The mean coverage or depth obtained ranged from 27 to 105 and the similarity base to base to the genome reference ranged from 98.5 to 99.6% (Supplementary Table S2). Both parameters indicate an acceptable quality to perform variant calling over the full genome.
Data deposition: the 87 genomes have been deposited at GenBank under the Accession Nos. MINA00000000.1, MTZQ00000000, MTZS00000000, MTZR00000000, and under submission SUB4397761 to be released upon publication (Bioproject Accession No. PRJNA485121) see also Supplementary Table S1.

SNP Identification and Selection
The comparison of the 89 genomic sequences led to the identification of SNP mutations. SNP were screened using information in the whole genome including intergenic, synonymous, and non-synonymous mutations.
SAMtools and BCFtools were used to transform the resulting BAM files format into VCF files. Then three tests were used to identify reliable SNPs; SNPs have a minimum of SAMtools quality score of 150. The number of forward (reverse) reads mapping with good quality onto the SNP site having the same base as in the reference had to be less than 20% of the total number of forward (reverse) reads mapping with good quality onto the SNP site having the variant base and the minimum forward (reverse) coverage on a SNP site has to be greater than 5. A total of 123 SNPs found in repetitive regions were filtered out. To perform this task, a window of 50 bp around the SNPs sites were blasted against the reference genome. The SNPs corresponding to sequences with more than one blast-hit (>99% match) were filtered out.
Single nucleotide polymorphisms distant on less than 10 bp from one another were excluded as they may be due to sequencing mistakes, as described before (Gardy et al., 2011).

Phylogenetic Analysis
Multilocus variable number of tandem repeat analysis and spoligotype patterns were carried out using Bionumerics 7.6.2 software (Applied Maths, Sint-Martens-Latem, Belgium) by drawing minimum spanning trees (MSTs) and dendrograms [Categorical (values) distance represented by UPGMA clustering method for MLVA or spoligotype, Multi-state categorical correlation represented by UPGMA clustering method for MLVA plus spoligotype]. SNP data were concatenated, resulting in a single character string (nucleotide sequence) for each strain. Phylogenetic and molecular evolutionary analyses were conducted using MEGA (version 7 2 ), using the neighborjoining method or maximum likelihood method with 1,000 bootstrap replicates, with distance calculated using the number of different SNPs.

SNP on Virulence Genes
We investigated whether it was possible to identify virulence traits related to the different phylogenetic lineages of French M. bovis by SNP analysis. To do this, we searched for SNP mutations among 181 genes that were described to be essential for virulence of the MTBC. The list of genes included in this study (Supplementary Table S4) were selected from the studies of Sassetti and Rubin (2003); Joshi et al. (2006), and Forrellad et al. (2013). Based on their function, molecular features or cellular localization, these 181 target genes included the following categories (Supplementary Table S4): Lipid and fatty acid metabolism (n = 25), Catabolism of cholesterol (n = 8), Cell envelope proteins (n = 12), 'mce' families proteins (n = 21), Lipoproproteins (n = 6), Metal-transporter proteins (n = 9), PE and PE_PGRS families regulators (n = 12), Secretion system (n = 17), Proteases (n = 6), Genes and metabolisms regulation expression (n = 26), Defense mechanisms (n = 30), and Other unknown function (n = 9).
The sequence reads of each strain were mapped to the assembled genome of M. bovis AF2122/97 to determine a similarity coefficient based on categorical difference. The clustering method UPGMA was used to depict the evolutionary relationships between the French M. bovis field isolates' genomes and the other genomes of the MTBC members presented in Figure 1.
Consistent with the predicted evolutionary scenario, the phylogenetic analysis showed M. canettii, is the most distantly FIGURE 1 | Whole genome-based evolutionary analysis of the panel of French field isolates of Mycobacterium bovis. The phylogenetic tree with numbered above the branches indicating the similarity coefficient between the organism and its common ancestor was constructed with whole genome, using the software bionumerics vs. 7.6.2 with the settings: similarity coefficient by categorical, scaling factor of 100 and clustering method by UPGMA.
related to the M. bovis strains with a coefficient of 98 followed by M. tuberculosis and M. africanum, the two preferentially human strains. M. bovis field isolates are closely related, nonetheless well-differentiated groups were distinguished. The BCG strains and those of the SB0120 group that share the same spoligotype appeared to be closely related and more distant to the strain AF2122/97 used as reference.

WGS and SNP Analysis
Whole-genome sequences were analyzed in order to detect SNP in comparison with the reference genome of M. bovis AF2122/97 allowing us to identify possible phylogenetic lineage markers and thus to resolve the phylogeny within M. bovis field isolates. The new 87 genomes sequences of the M. bovis field isolates plus the sequences of reference strains (M. bovis Pasteur_1173P2 BCG; M. bovis Tokyo_172 BCG) were mapped to the M. bovis reference strain (AF2122/97-BX24833) to identify SNPs.
A total of 9,170 SNPs were identified, 87% were located in coding regions, 2,880 of them were silent and 5,118 were missense. The distribution of SNPs covers the whole chromosome as depicted in Figure 2 and detailed in Supplementary Table S3.
In previous works, six groups, including the Eu 1 and Eu 2 clonal complexes, were distinguished within the population of M. bovis strains present in France (Rodriguez-Campos et al., 2012Smith, 2012;Hauer et al., 2015). These groups were defined by the combination of their spoligotypes, MLVA patterns and other mutations (summarized in Table 2).
The established phylogenetic analysis based on SNPs categorized all M. bovis isolates into groups which are congruent with those previously established with the combination of spoligotyping and MLVA (Figures 3A,B).
The results presented by the MST in Figures 3C,D show that the group resolution is better with SNPs analyses in comparison to combination of spoligotyping and MLVA. These observations suggest that SNP are better phylogenetic markers for M. bovis.
Not surprisingly, we observe that the highest homoplasy rate is given by spoligotyping (0.647) followed by MLVA (0.578) and significantly lower (0.004) by the use of SNPs (Supplementary Figure S2).
In addition, these analysis show that all the strains which clustered in the European 2 complex possessed the mutation in guaA (G→A; genome position 3765573 and gene position 843) unique of this complex (Rodriguez-Campos et al., 2012).    Table 1 are represented by name and their associated color.
were identified, in at least one strains, in our study. The details of their position and the genes impacted are indicated in the Supplementary Table S8.

Phylogenetically Informative SNPs
The dendrogram presented in Figure 4 details the numbers of SNP separating each sequenced strain distributed in the nine clusters, six of which have been previously described in the literature (Rodriguez-Campos et al., 2012;Smith, 2012;Hauer et al., 2015) and three less representative new groups. On a total of 9,170 SNPs, we established a catalog of SNPs specific to each clonal family (Figure 4 and Supplementary Table S3).
Cluster A (F4 group), previously defined by the lack of spacer 33 in the DR region and by a truncated repetition of the QUB26 locus, possesses 74 specific SNPs and is closely related to a new cluster B, defined by 52 specific SNPs, which is composed by only three strains from our selection. The next group on the phylogenetic tree is C (SB0134 group), which was previously characterized by the absence of the spacer 4 and 5 in the DR region, is only defined by seven SNPs but seems to be divided in two well-defined (respectively 28 and 29 specific SNPs) and balanced (respectively 4 and 5 strains) sub-groups. The last cluster on this branch includes the Eu1 complex, previously characterized by the absence of the spacer 11 and the RDEu1 deletion, which is defined here by 58 specific SNPs and composed by six of our strains and the AF2122/97 reference sequence. Another branch of the phylogenetic tree is composed by two cluster: a new one, "E, " only composed by three strains and defined by only six specific SNPs, and the well-defined (95 specific SNPs) Eu2 complex, also characterized by the absence of the spacer 21 in the DR locus and the guaA mutation (genome position 3765573 and gene position 843). The last and wellseparated branch (119 SNPs) is subdivided into three different clusters. Cluster G (F9 group), defined by 98 SNPs, is closely related to a new cluster, cluster H, which is composed by only two strains from our selection. Cluster I (SB0120 group) is only separated from the two previous ones by five SNPs. This cluster is composed by 30 field strains and two BCG reference sequences. The reference sequences clustered with one field strain (A7) and are well-separated from the other field strains, which form a clear group defined by 66 SNPs. This analysis confirms the coexistence of several groups in France, some of them already defined as clonal groups. Three of them represent almost 90% of the French strains: group I (SB0120), the largest one with 44.7% of the strains, and group A (F4 family) and group C (SB0134 group), with 28.5% and 14.5% respectively. Three other groups had been described before but are less preponderant in the French population: the Eu2 complex, dominant in the Iberian Peninsula, represents 9% of the French strains; the Eu1 complex, dominant in the British Islands, represents 1% of the French diversity; finally the group G (F9 family), a small group representing only 0.3% of the French diversity but clearly defined by a particular spoligotype pattern (absence of spacer 1-17 and 23-34) and localized in the south west of France. The new groups B, E, and H are not defined by any other characteristics on their spoligotype pattern or MLVA markers but regarding the phylogenetic analysis they seem to be clearly separated from the other groups. Further analysis on a larger set of strains is necessary to better define these groups. These results suggest that SNPs that are able to distinguish each cluster can be used to develop a novel SNP-based MLST to facilitate global epidemiology studies.

SNP: From Phylogeny to Virulence Traits
In order to learn if the above described groups could present differential virulence traits, the non-synonymous SNP (nsSNP) variants were investigated on 181 virulence genes (Supplementary Table S5). On these 181 genes, 30 were not affected in any strain genome (Supplementary Table S6), 151 genes have at least one SNP and in total 568 SNP, 176 synonymous and 392 non-synonymous, were identified (Supplementary Table S5). An overview of the distribution of the non-synonymous SNP (nsSNP) detected in each strain, organized by group, is illustrated in Figure 5. Three situations could be observed. In the first one, SNPs are distributed randomly (most frequent case see Figure 5). In the second situation, the same mutations are found in all French isolates at the same position (n = 4) (four full rays in the Figures 5, 6B) indicating that these four SNP are specific of the AF2122/97 reference strain. The third case, which is the most interesting one, shows mutations that are specific to either several or only just one family (partial concentric line). As shown in Figure 6B the vast majority of nsSNPs were detected in only one genome (236/392). Some genes represent a higher degree of mutation (rate median = 2) ( Figure 6A). This is particularly the case of six genes (Supplementary Table S7): Mb2074c (pKS12, 12457 amino acid) with 49 SNP variants, Mb3043c (PPE46, 1315 amino acids) 22 variants, Mb1849c (PE-PGRS33, 1507 amino acids) 14 variants, Mb2418 (PE_PGRS41, 1084 amino acids) 13 variants, Mb 1679c (PE_PGRS30, 3058 amino acids) nine variants and Mb1689 (pks7, 6382 amino acids) eight variants. Except for the Mb2074c, the rate of detected mutations is not always correlated to the size of the gene.
The dendrogram shown in Figure 7 highlights that nsSNPs analysis targeting virulence genes is also able to discriminate strains by group. Collectively, these results suggest that in addition to their phylogenetic differences, each group may have differential phenotypic traits.

DISCUSSION
In our previous study, we observed that the evolution of M. bovis strains responsible for tuberculosis outbreaks in France was accompanied by a significant reduction in their genetic diversity (Hauer et al., 2015). Indeed, current bTB cases in both livestock and wildlife are due to the spread of three predominant local genotypes. These observations were established by using the classical genotyping methods of spoligotyping and MLVA based on MIRU-VNTR (Hauer et al., 2016). While these targets and tools are usually sufficient for molecular epidemiology studies, they represent small hypervariable regions within the genome that generally evolve at a higher rate than the rest of the genome. Thus, in this study, we investigated the use of SNPs to define the genetic diversity extent in M. bovis that may provide insights into the evolution, transmission, molecular epidemiology and pathogenicity of the M. bovis isolates previously defined.
We selected 87 field isolates of M. bovis from pre-established clonal complexes. In order to identify informative SNPs among these genomes, whole genomes were sequenced and compared, while mutations were exhaustively searched.
Although these strains are genetically very close, more than 9,000 SNPs have been identified by genome comparison with commonly used filters. In some studies, SNPs were identified from a selection of genes of known functions (Behr et al., 2000). In the current study, SNP-based phylogenetic analysis using the whole set of SNP is congruent with the clusters defined by spoligotyping and MLVA.  Even if phylogeny results of genotyping (combination of spoligotyping and MLVA) and SNP based analyses are concordant, some differences may exist. For instance, strains belonging to the F9-family are closely related to the SB0120 group while conventional spologtyping-MLVA linked them to the Eu1 complex. Strains previously described as part of  F4-family by spoligotype and MLVA were confirmed by SNPs analyses. A previously defined complex, European 2, is characterized by a special mutation that appears mostly in the same SNPs group.
However, to test the performance and the degree of resolution of these SNP-based typing, the panel of strains must be adapted. In this study, strains were chosen from our collection from 1983 to 2011 to represent as closely as possible the genetic diversity of French M. bovis strains. With other different set of strains, SNPbased genotyping should be able to resolve misclassification of clonal clusters, to identify patterns of host or spatial association, and to specify their phylogeny and genealogy.
Our results strongly suggest that many new SNP targets can be identified by comparing the genomes of field M. bovis strains and may be used for the development of high-resolution genotyping methods as well as other subjects in relationship to strain evolution and their transmission. This set of about 9,000 SNPs has to be reduced and detailed for a wider use.
Non-synonymous SNPs targeting virulence genes were studied to reveal potential specific phenotypic traits within wellseparated phylogenetically French M. bovis isolates. Our analysis made it possible to identify non-synonymous SNPs in virulence genes unique to certain groups.
The two other genes presenting most mutations (Mb1689 pks7 and Mb2074c pks12) encode polyketide synthase. These modular enzymes are required for the production of unique lipids or glycolipid conjugates, which are critical for virulence, and also for other specific components of the complex mycobacterial cell envelope. It has been reported that pks7 plays an important role in virulence during the persistent phase of infection (Rousseau et al., 2003). Pks12 was shown to encode two sets of domains needed to produce fatty acids involved in phthiocerol synthesis, the diol required for dimycocerosyl phthiocerol (DIM) synthesis (Rousseau et al., 2003;Matsunaga and Sugita, 2012).
The non-synonymous SNP detected in virulence genes are most likely to impact on protein function and contribute to phenotypic variation. Further investigations on secreted proteins, epitope recognition, or pattern of outer membrane lipids (Yruela et al., 2016;Jankute et al., 2017;Ates et al., 2018) are needed to address the consequence of these mutations in order to define phenotypic differences across the main M. bovis groups still responsible for the persistence of outbreaks in cattle and wild life.
Sequential chromosomal nucleotide substitutions are the main driver of M. bovis genome evolution if we consider that the horizontal gene transfer seems not playing a role for the evolution of contemporary modern strain (Smith, 2012;Allen et al., 2013). Likewise, the genetic changes events that can be attributed to the transposable insertion sequences such as the IS6110 was not shown to affect the core genome of M. bovis (Otal et al., 2008;Allen, 2017).
Differences in genome sequence, including gene families that are important for bacterial infection and transmission, highlight differences with putative functional implication between isolates otherwise classified within the same spoligotype group.
A means to refine the clonal complex identity of the different groups identified in this study could include the analysis on the strains' pangenomes to detect the presence or absence of large deletions or insertions. To do this, it is necessary to create new reference genome sequences obtained by de novo assemblies of combination of second and third generation sequencing of strains belonging to main clonal groups. This would enable us to perform studies which are not limited to the core genome and to improve our knowledge on specific phenotype traits of M. bovis strain, the dynamics of strain transmission and evolution in multi-host systems.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

FUNDING
This project was funded by the French General Directorate for Food (DGAL) through the Réseau Français de Santé Animale.

ACKNOWLEDGMENTS
We thank the French General Directorate for Food (DGAl) for providing us with bTB epidemiological data. We are grateful to Victoria Boschiroli for her useful comments on the article.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2019.00955/full#supplementary-material FIGURE S1 | Geographical distribution of the 87 sequenced field isolates of Mycobacterium bovis. The strains were selected on the basis of their spoligotype and MLVA profile that, in different colors, represent all the genetic diversity of French M. bovis strains. This panel includes strains isolated from cattle (circles), goats (triangles), sheep (stars) and wildlife (losanges), pigs (pentagons), wild boar (squares).
FIGURE S2 | Comparison of the homoplasy index among the different genotyping methods. HI was calculated with software package PAUP based on the number of observed changes at each character compared to the expected number of changes assuming absence of homoplasy.