Eight-barbel loaches belonging to the genus Lefua (Balitoridae, Cypriniformes) are primary freshwater fish less than 10 cm in total length. The evolutionary process in Lefua forms part of the evolutionary history of Japanese and East Asian freshwater fishes. Four species have been formally described. Lefua nikkonis (Jordan et Fowler, 1903) is endemic to Hokkaido, Japan, L. echigonia Jordan et Richardson, 1907 inhabits Tohoku to Kinki districts in Honshu, Japan, and L. costata (Kessler, 1876) is distributed widely in East Asia including Korea, China, and Russia, but has been introduced to restricted areas in Honshu, Japan  . Lefua pleskei (Herzenstein, 1887) is distributed in eastern Russia  . Another species, Lefua sp. was separated from L. echigonia based on morphological traits with a greater distance between the dorsal and ventral fins, longer snout length, lower body height, and narrower body width of Lefua sp. than L. echigonia   , and on ecological traits with habitats of relatively fast-flowing mountain streams with gravelly beds of Lefua sp. contrasting with habitats of relatively slow-flowing streams with muddy beds in marshlands, spring water, rice paddy irrigation channels, and the backwaters of floodplains of L. echigonia  . Lefua sp. is waiting for a formal scientific description. This species is distributed in western Japan including Kinki, Hokuriku, and Chugoku districts in Honshu and also in Shikoku. Additionally, we suggested recently that Lefua loach occurring locally in Aichi and Shizuoka prefectures of the Tokai region in Honshu, Japan is the sixth possible species   . Although the loach had been identified as Lefua sp. on morphological and ecological grounds as above, our phylogenetic study by mitochondrial DNA (mtDNA) and nuclear DNA (nDNA) sequencing showed that the loach in the Tokai region comprised a monophyletic group and was more closely related to L. echigonia than to Lefua sp. In our previous study, we designated it the L. sp. Tokai population, and thereafter Lefua sp. in western Japan was classified as Lefua sp. 1 and the loach in the Tokai region as Lefua sp. 2 in the Red List issued by the Japanese Ministry of the Environment in 2014 (see also  ). According to the classification, we herein designate the loach in the Tokai region as Lefua sp. 2. Some studies showed embryological differences in Lefua loaches     . However, we cannot perform comparative embryology among the above six species and among intraspecific populations of each species, because those studies have been done using only limited samples of two species, Lefua sp. 1 and L. echigonia.
Although the taxonomic status of Lefua sp. 2 remains to be specified, it has a unique evolutionary background clearly  . We suggested parallel morphological evolution between Lefua sp. 1 and Lefua sp. 2  . Although their distributions are disjunct, both species occur in relatively fast-flowing mountain streams with gravelly beds and generally have more slender bodies for adapting to fast- flowing streams than L. echigonia. The slender bodies allow them to avoid fast flowing water through exploitation of spaces within the gravel bed of streams. Lefua sp. 2 and L. echigonia inhabit the Tokai region. Complex geological events in their distribution range are likely associated with the evolutionary process. Upliftment in the whole range might have forced ancestors of Lefua sp. 2 to adapt to fast-flowing mountain brooks  . Later L. echigonia invaded the Tokai region during erosion of the western flanking zone (Ryoke metamorphic belt) of the Median Tectonic Line (MTL). At present, Lefua sp. 2 preferentially inhabits the eastern flanking zone (Sanbagawa metamorphic belt) of MLT and the zone of uplifted volcanic and marine deposit. The Ryoke metamorphic belt tends to be strongly eroded, producing relatively flat terrain suitable for development of L. echigonia habitat. This is not the case in the Sanbagawa metamorphic belt of a different rock composition.
Because of its extinction risk and unique evolutionary background, Lefua sp. 2 was registered as endangered species in the Red List as well as L. echigonia, Lefua sp. 1, and L. nikkonis    . These Lefua loaches have been severely threatened by habitat destruction due to human activities. We demonstrated that L. echigonia comprised genetically seven intraspecific populations (Hokuriku, Tohoku, Yamagata, North-Kanto, South-Kanto, Kinki, and Tokai) and Lefua sp. 1 comprised two intraspecific populations (Kii-Shikoku and Sanyo)   . These populations occur in regions that are geographically well separated by mountain ranges and highlands   . The populations should be considered evolutionary significant units (ESUs, as defined by Moritz  ) for the protection of the endangered loaches  . However, intraspecific phylogenetic relationships and variation in Lefua sp. 2 have not been investigated thoroughly. Elucidation of cryptic ESUs is crucial for the protection of this unique loach.
In the present study, we demonstrate intraspecific phylogenetic relationships and variation in Lefua sp. 2 and L. echigonia (Tokai population) by sequencing the mitochondrial D-loop region. We present fundamental information for the protection of these endangered species, and discuss their evolutionary history in the Tokai region.
2. Materials and Methods
Collecting localities of Lefua sp. 2 (83 specimens) and L. echigonia (40 specimens) are shown in Figure 1 and information on all samples including outgroup specimens is listed in Table 1. Since both species was registered as endangered species in the Red List, we used large circles on Figure 1 to conceal detailed information of habitats for the conservation of the two endangered species. Specimens of Lefua sp. 1 obtained from Hidaka (Sanyo population) and L. echigonia obtained from Ouchi (Hokuriku population) were used as the outgroup for the Lefua sp. 2 phylogeny. Those of L. echigonia obtained from Aogaki (Kinki population) and Ouchi were used as the outgroup for the L. echigonia phylogeny in the Tokai region. Lefua sp. 2 and L. echigonia are sympatric in Sanagawa, Aichi Prefecture (sample Nos. 69 and 111), and occur in the close vicinity in Kosai, Shizuoka Prefecture (sample Nos. 80 and 118).
Table 1. Sample list.
*Deposited previously; **Deposited previously and revised in this study.
Figure 1. Locations in the Tokai region where specimens of Lefua sp. 2 and L. echigonia were collected Refer to Table 1 for details of the collection sites and sample numbers. We used large white circles for sample localities of Lefua sp. 2 and large black circles for those of L. echigonia. As a result, each circle does not always represent a single locality and in some cases integrated pleural localities. The numbers accompanied with circles denote the sample numbers in Table 1. A to I accompanied with circles include the following samples; A, 1-20; B, 22, 25-44; C, 45-50, 55, 57, 63; D, 51, 56, 58, 61, 66, 68; E, 69, 109, 110, 111; F, 60, 62, 64; G, 80, 82; H, 77, 81, 83; I, 91, 92.
2.2. DNA Sequencing
Total DNA from fin or muscle of each fish was prepared using a DNeasy® Blood & Tissue Kit (QIAGEN GmbH, Hilden) according to the manufacturer’s protocol. To amplify mtDNA containing the D-loop region, PCR was performed in a reaction solution (50 µl) containing template DNA (2 µl) and KOD dash (Toyobo Co., Ltd., Osaka) with an initial denaturation at 94˚C for 2 min, followed by 30 cycles of denaturation at 94˚C for 30 sec, annealing at 55˚C for 10 sec, and extension at 74˚C for 30 sec, and with final extension at 74˚C for 7 min. Sequences of the primers used for amplification are shown in Table 2  . The amplified DNA fragment was purified using a QIAquick Purification Kit (QIAGEN GmbH, Hilden). Direct sequencing of the double-stranded PCR product was performed using an ABI PRISM BigDye® Terminator v1.1 Cycle Sequencing Kit (Applied Biosystems Inc., California) and the primers used for PCR on Model 377 and 377XL DNA sequencers (Applied Biosystems) according to the manufacturer’s directions. Alternatively, direct sequencing was performed using a Genome Lab™ DTCS-Quick Start Kit on a CEQ™ 2000XL DNA Analysis system (Beckman Coulter Inc., California) according to the manufacturer’s directions.
New sequence data were deposited in the DDBJ database under accession numbers LC062792-062892. We had previously determined mtDNA sequences used in this study (AB102843-102846, AB102849, AB102850, AB177672, AB177677, AB177699, AB251866, AB251870-251878, AB599746, AB599766- 599770     ).
Table 2. Primers used in this study.
2.3. Phylogenetic Analysis
DNA sequences of the mitochondrial D-loop region were edited and aligned using DNASIS (Hitachi Software Engineering Co., Ltd., Tokyo) and MEGA 6.0  , and the alignments were corrected by visual inspection. We used 825 bp for the Lefua sp. 2 phylogeny and 867 bp for the L. echigonia phylogeny in the Tokai region, excluding indels and ambiguous sites. Neighbor-joining (NJ) and maximum parsimony (MP) trees were constructed using MEGA 6.0 and PAUP*4.0 beta10  , respectively. Genetic distances were computed by Kimura’s two- parameter model  . Tree reliability was evaluated by generating 1,000 bootstrap replicates. The majority-rule consensus MP tree was constructed by conducting a heuristic search based on the 1,000 bootstrap replicates with an unweighted ts/tv ratio. The Bayesian (BA) tree was constructed using MrBayes version 3.1.2  based on the model evaluated by the MrModeltest 2.3  . The best models were GTR + I + G for the Lefua sp. 2 phylogeny and HKY + G for the L. echigonia phylogeny. The Monte Carlo Markov chain (MCMC) length was 5 × 106 generations, and we sampled the chain every 100 generations. MCMC convergence was assessed by calculating the potential scale reduction factor, and the first 1 × 104 generations were discarded. The minimum spanning tree representing relationships among haplotypes was constructed using MEGA 6.0 and haplotype and nucleotide diversities and Genetic differentiation (Fst) were calculated using Arlequin 188.8.131.52  . The statistical significance of Fst was evaluated by calculating 1 × 104 values.
1) Phylogenic relationships in Lefua sp. 2 and L. echigonia from the Tokai region
To determine phylogenetic relationships of 83 specimens in Lefua sp. 2, the NJ tree was constructed based on 825-bp sequences using Lefua sp. 1 from Hidaka and L. echigonia from Ouchi as the outgroup (Figure 2). There were 58 variable sites and 29 parsimony informative sites. Most branches were very short, indicating high sequence similarities. The monophyly of Lefua sp. 2 was well supported (NJ, 100; MP, 100; BA, 1.00), but other groupings within Lefua sp. 2 were generally not well supported because of high sequence similarities. We arbitrarily designated assemblies including specimens mainly from Yahagi, Toyo, and Tenryu River systems as Groups 1 to 3, respectively. Group 1 comprised specimens collected in the major Yahagi and neighboring minor Mito, Nishida, and
Figure 2. Phylogenetic relationships in Lefua sp. 2 based on the 825-bp mitochondrial D-loop sequences. The NJ tree was constructed based on genetic distances calculated with Kimura’s two-parameter model using a total of 85 samples including Lefua sp. 1 from Hidaka and L. echigonia from Ouchi as the outgroup. The scale bar indicates 0.01 substitutions per site. Statistical supports of tree topologies are specified with NJ (left) and MP (middle) bootstrap values and BA posterior probabilities (right) in the vicinity of corresponding nodes. Simbols denote river systems: ▲, Yahagi; ∇, Mito; ?, Nishida; △, Otowa; ○, Toyo; ●, Sana; ■, Tenryu; □, Miyakoda; ◆, Oota.
Otowa River systems (hereafter abbreviated as the Yahagi River system), and was marginally supported (NJ, 69; MP, 49; BA, 1.00). Group 2 was consisted of specimens collected in the major Toyo and neighboring minor Sana River systems (hereafter abbreviated as the Toyo River system), and was paraphyletic to Group 3. Group 3 comprised specimens collected in the major Tenryu and neighboring minor Miyakoda and Oota River systems (hereafter abbreviated as the Tenryu River system), and was poorly supported (NJ, 3; MP, 14; BA, -). Group 2 + Group 3 is monophyletic, albeit poorly supported. Group 2 included exceptionally five specimens collected in the Tenryu River system. As described in the Discussion, localities of those specimens from the Tenryu River system, Gotengawa, Takihashigawa, Fukayagawa, Nanegawa, and Nagaishigawa (sample Nos. 70-74 in Table 1), are very close to tributaries of the Toyo River system. Group 3 included an exceptional specimen from the Toyo River system (Aderagawa, sample No. 65).
The p-distance and Fst were low between Groups 2 and 3, although those between Groups 1 and 2 and between Groups 1 and 3 were relatively higher (Table 3), suggesting that Lefua sp. 2 in the Yahagi River system was moderately differentiated from that in the other river systems. MP and BA trees demonstrated similar basal divergences, but splitting at the tip was not conservatively recovered (data not shown).
To determine phylogenetic relationships of 40 specimens in L. echigonia, the NJ tree was constructed based on 867-bp sequences using L. echigonia from Aogaki and Ouchi as the outgroup (Figure 3). There were 105 variable sites and 26 parsimony informative sites. The monophyly of specimens from the Tokai region, i.e. the Tokai population of L. echigonia, was well supported (NJ, 95; MP, 97; BA, 0.99). We arbitrarily designated assemblies as Groups A and B. Group A comprised specimens from Shizuoka and southeastern Aichi prefectures (NJ, 41; MP, -; BA, 0.90) and Group B from northwestern Aichi, Gifu, and Mie prefectures (NJ, 47; MP, 56; BA, 0.61). The p-distance and Fst between Groups A and B were low, but their p-distance was higher than those within each group (Table 4). MP and BA trees demonstrated similar basal divergences, but splitting at the tip was not conservatively recovered (data not shown).
2) Haplotype relationships in Lefua sp. 2 and L. echigonia from the Tokai region
Forty-five haplotypes were detected in Lefua sp. 2, and the minimum spanning tree was constructed to represent their relationships (Figure 4). The tree showed a complicated network of haplotypes, but haplotypes represented by
Table 3. Genetic divergence in Lefua sp. 2.
Above diagonal, Fst; on diagonal, intragroup p-distance; below diagonal, intergroup p-distance. *statisti- cally significant (P < 0.05).
Figure 3. Phylogenetic relationships in L. echigonia based on the 867-bp mitochondrial D-loop sequences The NJ tree was constructed based on genetic distances calculated with Kimura’s two-parameter model using a total of 42 samples including L. echigonia from Aogaki and Ouchi as the outgroup. The scale bar indicates 0.005 substitutions per site. Statistical supports of tree topologies are specified with NJ (left) and MP (middle) bootstrap values and BA posterior probabilities (right) in the vicinity of corresponding nodes.
Figure 4. Minimum spanning tree of 45 haplotypes in Lefua sp. 2 Circles indicate haplotypes and their dimensions are proportional to the number of individuals possessing the same haplotype. Black circles on lines indicate hypothetical haplotypes not found in this study. The numbers denote the sample numbers in Table 1. Dark grey (purple in the color version), Group 3; pale grey (yellow), Group 2; intermediate grey (green), Group 1.
Table 4. Genetic divergence in L. echigonia.
Above diagonal, Fst; on diagonal, intragroup p-distance; below diagonal, intergroup p-distance; *, statistically significant (P < 0.05).
specimens belonging to the three groups were roughly separated from one another. The haplotypes of those belonging to the Yahagi River system were relatively well separated from the others, supporting that Lefua sp. 2 in the Yahagi River system was genetically differentiated from the others. The haplotypes of Gotengawa, Takihashigawa, Fukayagawa, Nanegawa, and Nagaishigawa from the Tenryu River system were included in the assembly of haplotypes from the Toyo River system, while those of Aderagawa, Fumagawa (sample No. 81), and Morimachiichimiya (83) from Group 3 were included in the assembly of haplotypes from Group 2. Therefore, Lefua sp. 2 in the Toyo River system was not genetically distinguishable well from that in the Tenryu River system.
Twenty-eight haplotypes were detected in L. echigonia, and the minimum spanning tree was constructed to represent their relationships (Figure 5). A star-like structure of haplotype relationships was represented with the haplotype of the greatest majority shared by 10 specimens from eastern Aichi prefecture.
3) Haplotype and nucleotide diversities in Lefua sp. 2 and L. echigonia from the Tokai region
Values of the haplotype diversity in Lefua sp.2 were relatively high and similar between groups (Table 5), but values of the nucleotide diversity were very low and those in Groups 1 and 2 (Yahagi and Toyo River systems) were lower than that in Group 3 (Tenryu River system). High haplotype and low nucleotide diversities in L. echigonia were similar between Groups A and B.
Our results showed that Lefua sp. 2 was an established entity in the Lefua loaches (see also   ), and that there were not well-defined genetic structures discernible within Lefua sp. 2. Taking the D-loop as the fastest evolving region in mitochondrial DNA into account, Lefua sp. 2 dispersed quite recently, or it has been maintained as a single population for long time. Lefua sp. 2 is more closely related to L. echigonia than to Lefua sp. 1. Divergence of Lefua sp. 2 and L. echigonia was estimated at 1.4 - 1.5 million years ago, whereas that of the Tokai and Kinki populations of L. echigonia at 0.7 million years ago   . Despite its long history, Lefua sp. 2 represented low intraspecific genetic distances (Table 3) and low nucleotide diversity, although haplotype diversity was relatively high because a number of haplotypes with small nucleotide changes were detected (Table 5). Therefore, Lefua sp. 2 might have once reduced its population size and dispersed recently from the relic area(s), which is not the case when its
Figure 5. Minimum spanning tree of 28 haplotypes in L. echigonia Circles indicate haplotypes and their dimensions are proportional to the number of individuals possessing the same haplotype. Black circles on lines indicate hypothetical haplotypes not found in this study. The numbers denote the sample numbers in Table 1. Dark grey (pink in the color version), Group A; pale grey (blue), Group B.
Table 5. Haplotype and nucleotide diversities in Lefua sp. 2 and L. echigonia.
n, number of samples.
evolutionary rate has been extremely slow. However, specimens from the Yahagi River system were moderately differentiated from those from the Toyo and Tenryu River systems, although it is difficult to distinguish definitely specimens between the Toyo and Tenryu River systems. Thus, we can consider that Lefua sp. 2 consists of two natural populations, that is, the Yahagi River population and the Toyo-Tenryu River population.
Star-like relationships of haplotypes in L. echigonia suggested the dispersal origin located in eastern Aichi prefecture. This does not look consistent with our assumption that L. echigonia invaded the Tokai region from the western habitat, near the boundary between the distributional ranges of the Tokai and Kinki populations. The invasion was caused by erosion of the western flanking zone of MLT and the development of relatively flat terrain more suitable for L. echigonia than for Lefua sp. 2. However, it is possible that the dispersal origin of the former western habitat has already been lost, possibly because of destruction by human activities.
In Lefua sp. 2, the specimens from Gotengawa, Takihashigawa, Fukayagawa, Nanegawa, and Nagaishigawa in the northwestern Tenryu River system were more closely related to those in the Toyo River system than to those in the southeastern Tenryu River system. The distance between tributaries of the northwestern Tenryu River and northeastern Toyo River is less than 1.5 km. Therefore, we assumed that Lefua sp. 2 was introduced from the Toyo River to Tenryu River by river capturing. This provides a biological clue in resolving whether river capturing between the Tenryu and Toyo Rivers occurred previously. River capturing has been suggested based on river trajectories, while Ikeda  insisted from research of riverine sediments that there was no river capturing between the two rivers. The present study indicated that at least a northeastern part of the Toyo River was captured by a northwestern part of the Tenryu River. The specimen from Aderagawa in the Toyo River system was exceptionally included in Group 3. The collection site of Aderagawa is very close to the tributary of Tenryu River. This also suggests river capturing between the Toyo and Tenryu Rivers, although we cannot eliminate the possibility of artificial transplantation of Lefua sp. 2 from the Tenryu River to Toyo River. Nevertheless, the present results suggest that the Toyo River and Tenryu River systems form a single population through gene flow facilitated by river capturing. Population genetic studies of various taxa can present more solid information of geological settings such as river capturing in the Tokai region.
The two populations should be considered as evolutionary significant units (ESUs) for the conservation of Lefua sp. 2. The translocation of fishes is severely prohibited by a guideline issued by the Ichthyological Society of Japan; however, there are some problems to follow this guideline. One of them is how to define ESUs. We consider it important to regard a natural population as an ESU. Many fish habitats have been fragmented by anthropogenic activities. Especially for fishes inhabiting montane streams, such as Lefua sp. 2, dam and weirs construction is critical by preventing fishes from returning upstream to their original habitats when they are displaced downstream by strong currents after heavy rains. Individuals are often isolated from each other and divided into small assemblages, and thus a natural population is destroyed by the obstruction of gene flow. Small assemblages are susceptible to bottlenecking, and haplotypes included in an original population can be fixed in different manners in different tributaries. When assemblages are investigated using genetic makers, they can be recognized as genetically distinct and thus reasoned to be assigned to different ESUs. However, they are small anthropogenic populations, but not natural populations. In this case, translocation is prohibited based on anthropogenic but not natural grounds. Next is how finely ESUs should be defined. We still have genetic markers separating organisms up to the individual level, and thus we can define ESUs more finely using faster-evolving genetic markers. However, it is very difficult or impossible and even impractical to conserve so many finely defined ESUs. We have to resolve these problems to adopt the best procedure for the fish conservation. Specimens of L. echigonia from the Tokai region showed lower genetic differences between Group A and B (0.008 in Table 4) than those of Lefua sp. 2 between Group 1 and Group 2 or 3 (over 0.010 in Table 3). Presently we cannot say whether two groups of L. echigonia can be regarded as ESUs.
Many Japanese freshwater fishes are facing to extinction. Urgent measures appropriate to each fish species are indispensable to protect them and guide lines are needed to be further refined. The present study contributes to protection of endangered eight-barbel loaches by giving fundamental genetic information and defining new ESUs in Lefua sp. 2. Studies such as the present study are crucial and will together lead to comprehensive conservation of Japanese freshwater fishes in future.
We express our sincere thanks to Mr. Michikazu Matsuda, Mr. Masayoshi Sugiura, Dr Seiichi Mori, Dr. Kazumi Hosoya, Mr. Tesuo Kamei, Mrs. Yumiko Yamashina, Mrs. Keiko Muraoka, Mr. Osamu Inaba, Mr. Koichi Kirihara, and Mr. Sakae Nakamura for their assistance in collecting loaches and for providing useful information on the biology of Lefua loaches.
 Mihara, M., Sakai, T., Nakao, K., Martins, L. de O., Hosoya, K. and Miyazaki, J.-I. (2005) Phylogeography of Loaches of the Genus Lefua (Balitoridae, Cypriniformes) Inferred from Mitochondrial DNA Sequences. Zoological Science, 22, 157-168.
 Shedko, S.V., Miroshnichenko, I.L. and Nemkova, G.A. (2008) On the Systematics and Phylogeography of Eight-Barbel Loaches of the Genus Lefua (Cobitoidea: Nemacheilidae) mtDNA Typing of L. pleskei. Russian Journal of Genetics, 44, 817-825.
 Miyazaki, J.-I., Nakao, K., Mihara, M., Sakai, T., Gunji, Y., Tojo, K., Muraoka, K. and Hosoya, K. (2007) Incongruence between mtDNA Phylogeny and Morphological and Ecological Characters in Loaches of the Genus Lefua (Balitoridae, Cypriniformes). Zoological Science, 24, 666-675.
 Miyazaki, J.-I., Dobashi, M., Tamura, T., Beppu, S., Sakai, T., Mihara, M. and Hosoya, K. (2011) Parallel Evolution in Eight-Barbel Loaches of the Genus Lefua (Balitoridae, Cypriniformes) Revealed by Mitochondrial and Nuclear DNA Phylogenies. Molecular Phylogenetics and Evolution, 60, 416-427.
 Aoyama, S. and Doi, T. (2011) Morphological Comparison of Early Stages of Two Japanese Species of Eight-Barbel Loaches: Lefua echigonia and Lefua sp. (Nemacheilidae). Folia Zoologica, 60, 355-361.
 Sano, I., Kirikane, M., Ogata, K., Kobayshi, K., Tojo, K. and Miyazaki, J.-I. (2016) Phylogeny in Eight-Barbel Loaches of the Genus Lefua (Balitoridae, Cypriniformes) Based on Nucleotide Sequences of the Mitochondrial Cytochrome b Gene. Bulletin of the Faculty of Education (University Yamanashi), 1, 183-192. (In Japanese)
 Sakai, T., Mihara, M., Shitara, H., Yonekawa, H., Hosoya, K. and Miyazaki, J.-I. (2003) Phylogenetic Relationships and Intraspecific Variations of Loaches of the Genus Lefua (Balitoridae, Cypriniformes). Zoological Science, 20, 501-514.
 Tamura, K., Stecher, G., Peterson, D., Filipski, A. and Kumar, S. (2013) MEGA6: Molecular Evolutionary Genetics Analysis Version 6.0. Molecular Biology and Evolution, 30, 2725-2729.
 Kimura, M. (1980) A Simple Method for Estimating Evolutionary Rates of Base Substitutions through Comparative Studies of Nucleotide Sequences. Journal of Molecular Evolution, 16, 111-120.
 Huelsenbeck, J.P., Ronquist, F., Nielsen, R. and Bollback, J.P. (2003) Bayesian Inference of Phylogeny and Its Impact on Evolutionary Biology. Science, 294, 2310-2314.
 Posada, D. and Buckley, T.R. (2004) Model Selection and Model Averaging in Phylogenetics: Advantages of the AIC and Bayesian Approaches over Likelihood Ratio Tests. Systematic Biology, 53, 793-808.
 Excoffier, L. and Lischer, H.E.L. (2010) Arlequin Suite Ver 3.5: A New Series of Programs to Perform Population Genetics Analyses under Linux and Windows. Molecular Ecology Resources, 10, 564-567.