Phylogeography of Ara militaris (Military Macaw): Implications for Conservation

: The Military Macaw ( Ara militaris ) is an endangered bird species with disjunct geographic distribution across the Neotropics, consisting of three recognized subspecies: One in Mexico ( A. m. mexicanus ) and two in South America ( A. m. militaris and A. m. bolivianus ). However, due to the limited phenotypic differentiation between these allopatric taxa, their taxonomic status has been the subject of debate. In this study, we explored mitochondrial DNA ( mtDNA ) variability to determine the phylogeographical pattern through phylogenetic and ecological modelling analyses. We also aimed to describe the evolutionary relationships of twelve A. militaris populations. We identiﬁed 41 haplotypes in the 300 bp region of the Cytochrome b ( Cyt-b ) gene of the mtDNA and low nucleotide diversity. The observed phylogeographic structure suggests the existence of two clades: One composed of A. m. militaris and A. m. bolivianus and another consisting solely of A. m. mexicanus . The A. m. mexicanus clade further divides into two recognized subclades: Sierra Madre Oriental and northeastern portion of the Sierra Madre Occidental. Ecological analyses revealed that the niche similarity between these lineages was lower than expected by chance. Additionally, results from low cross-prediction tests indicated that the two lineages have inhabited different environmental spaces since the Late Pleistocene. This divergence may be associated with a steep ecological gradient and contemporary geographical barrier. Based on our results, we suggest that at least the A. m. mexicanus has a divergent evolutionary history; therefore, it should be considered as a different evolutionarily signiﬁcant and management unit. We recommend that future conservation strategies in Mexico incorporate effective protection measures, including habitat preservation and the reduction of illegal trade, to ensure the preservation of viable populations.


Introduction
Phylogeographical studies have held a central role in comprehending various aspects of Neotropical bird species, encompassing gene flow patterns, hybridization, historical range fragmentation and expansion, and speciation [1][2][3][4][5].These studies have employed several approaches to demonstrate how genetic and morphological divergence is profoundly shaped by the spatial structure of ecosystems, the presence of barriers (geographic and ecological) that impede gene flow, and the impacts of Pleistocene climatic changes [6].These factors, whether acting individually or in concert, can propel ecological adaptations to different environments, resulting in ecological niche specialization and lineage divergence [1].Despite numerous studies having provided insights into the biogeographical history of diverse taxonomic groups across the Neotropics, the evolutionary history of several psittacids species (Aves: Psittaciformes) remains partially understood [7,8].
The Neotropics harbor the highest diversity of psittacids worldwide, encompassing at least 165 acknowledged species within the region [9,10].Recent phylogenetic and phylogeographic analyses have established the tribe Arini as a well-supported monophyletic clade [3,8,[11][12][13][14].Moreover, they also suggest that Pleistocene climatic fluctuations played a significant role in the diversification of psittacids by driving genetic and phenotypic divergence within species and facilitating speciation [13,15,16].Molecular studies on endangered parrot species have further unveiled substantial population structure, attributed to historical events as well as more recent environmental changes [17,18].These findings are particularly important as they have led to the re-evaluation of the taxonomy and conservation status of isolated populations of many endangered species [18,19].They also highlight the significance of conducting similar studies on other psittacids, considering that 31% of the species in the Neotropics are classified as endangered taxa [20,21].
Large-bodied macaws are indeed among the most endangered parrot species globally [21,22].Unfortunately, they have experienced rates of extinction in recent times, as exemplified by the Cuban species (Ara tricolor) [23].The principal factors contributing to their decline encompass illegal hunting, habitat loss, fragmentation, and the effects of climate change [24][25][26].While substantial research has been directed towards the phylogeography and population genetics of species classified of least concern internationally and facing the risk of extinction according to local norms (Mexico and Central American countries), such as the Scarlet Macaw (Ara macao) [27], and critically endangered species like Spix's Macaw (Cyanopsitta spixii) [28], other Mesoamerican species, such as the Military Macaw (Ara militaris), have received comparatively diminished attention [29,30].
The Military Macaw (Ara militaris) boasts an extensive geographic range that stretches from Mexico to South America (Figure 1).Within this range, it consists of three subspecies with allopatric distribution [29]: A. m. militaris is found from the northwest of Venezuela to southeast Peru; A. m. bolivianus is distributed from northern Bolivia to northwest Argentina; and A. m. mexicanus is restricted to Mexico, spanning across the biogeographic provinces of Sierra Madre Occidental, Sierra Madre del Sur, and Sierra Madre Oriental [31].It is noteworthy that A. militaris is not present in Central America, as this region is occupied by its sister species, A. ambiguus [9].
Despite A. militaris and A. ambiguus being historically regarded as a single species complex [29], recent phylogenetic studies have found strong support for the monophyly of both species, thereby recognizing them as distinct evolutionary lineages [29,30].Moreover, these studies have revealed the presence of a genetic structure within A. militaris [29,30].However, these studies were based on samples obtained from museum specimens and did not cover the entire distribution range of the species, particularly in Mexico.In this study, we increased the number of populations of A. m. mexicanus, utilizing contemporary samples that encompass the complete distribution range in Mexico.Furthermore, we have explored whether niche divergence accompanies evolutionary differentiation within these taxa, as suggested by previous studies conducted in other regions of the Neotropics [32][33][34][35].This type of information can enable us to identify priority conservation units and promote effective planning for the protection efforts for the Military Macaw in the future [21].Considering that this species exhibits a clear population decline (with fewer than 10,000 individuals remaining in the wild), it is classified by the International Union for Conservation of Nature as threatened (vulnerable status), and the Convention on International Trade in Endangered Species of Wild Fauna and Flora classifies it as endangered [36,37].
the future [21].Considering that this species exhibits a clear population decline (with fewer than 10,000 individuals remaining in the wild), it is classified by the International Union for Conservation of Nature as threatened (vulnerable status), and the Convention on International Trade in Endangered Species of Wild Fauna and Flora classifies it as endangered [36,37].In this study, we conducted an analysis of the phylogeographic structure of A. militaris to infer evolutionary relationships and examine patterns of genetic diversity and structure among its populations.Additionally, we estimated the divergence time between the identified clades to investigate the potential influence of past climatic events on the taxon's diversification.Given that ecological differentiation may arise from geographic disruptions to the species' optimal ecological niche, we hypothesized that distinct phylogenetic lineages would exhibit discernible patterns of niche segregation [35].Based on these findings, we identified conservation units [36] and discussed their conservation implications and priorities for the Military Macaw in Mexico.In this study, we conducted an analysis of the phylogeographic structure of A. militaris to infer evolutionary relationships and examine patterns of genetic diversity and structure among its populations.Additionally, we estimated the divergence time between the identified clades to investigate the potential influence of past climatic events on the taxon's diversification.Given that ecological differentiation may arise from geographic disruptions to the species' optimal ecological niche, we hypothesized that distinct phylogenetic lineages would exhibit discernible patterns of niche segregation [35].Based on these findings, we identified conservation units [36] and discussed their conservation implications and priorities for the Military Macaw in Mexico.

Sampling
In this study, we focused our taxonomic sampling on the A. militaris complex.We analyzed the entire distribution range of the species, including the areas for its three recognized subspecific taxa: A. m. mexicanus, A. m. militaris, and A. m. bolivianus (Figure 1).For A. m. mexicanus, we collected feathers from 35 localities across its geographic distribution in Mexico (permit number, SGPA/DGVS/05767/22, Secretaría de Medio Ambiente y Re-cursos Naturales).The feathers were collected at the base of the trees in nesting, feeding, and resting sites during the breeding season.In some localities, feathers were directly collected from the nests (see Table 1 and Table S1).We considered only one feather per tree or per nest as a single individual per locality.Once collected, we cleaned the feathers with 90% ethanol and placed them in paper bags for transportation to the laboratory [37].Additionally, we obtained samples of the three subspecies of Military Macaw from 32 museum specimens housed in the American Museum of Natural History (Manhattan, NY, USA), by dissecting 20 to 30 mm 2 of toe pad skin (Table 1 and Table S1).Each of these samples was stored in 70% ethanol for transportation to the laboratory.In total, 67 samples from 23 locations were obtained.All sampling localities separated by less than 100 Km and within the same geographical region were considered a population, resulting in the recognition of 12 populations along the species distribution range (Table 2).Total genomic DNA was extracted from feather samples using a standard proteinase K/SDS digestion method, following the chloroform:ethanol purification protocol described by Leeton and Christidis [38].For skin samples, the Dneasy ® extraction kit (Qiagen, Hilden, Germany) was used, following the manufacturer's protocol.Given the challenges associated with working with old museum skin samples, which often yield small quantities of degraded DNA and are susceptible to contamination from foreign DNA, we followed the recommendations by Murphy et al. [39].These recommendations were implemented to enhance DNA quality and yield while minimizing the risk of contamination.
We conducted amplification and sequencing of fragments from the mitochondrial Cytochrome-b gene (Cyt-b), a commonly employed marker in phylogenetic and phylogeographic studies of parrots, e.g., [39][40][41].To amplify the fragments, we utilized the standard Cyt-b primers L14996 (5 -ATCTCHGCHTGATGAAAYTTYGG-3 ) and H15646 (5 -GGGGTGAAGTTTTCTGGGTC-3 ) designed by Sorenson et al. [42], resulting in fragments of 600 bp in length.However, for museum skin samples, we could only amplify shorter Cyt-b fragments (300 bp) using the internal primer L15413 (5 -GAATGAGCATGAG GCGGATTCTCA-3 ) in conjunction with H15646.The degradation of DNA in these museum samples, which were up to 124 years old, likely contributed to the limitation in fragment length.Nevertheless, the successful utilization of short mitochondrial (mtDNA) fragments has been demonstrated in other avian phylogeographic investigations, e.g., [39,[43][44][45], providing valuable and reliable information regarding population history [46].
For all PCR reactions, we employed the QIAGEN Multiplex PCR kit (Qiagen, Hilden, Germany).Each reaction comprised 2 µL of Master Mix, 0.5 µL of primers (5 pmol/µL), 1 µL of UV-sterilized distilled water, and 1 µL of total genomic DNA (20-50 ng/µL), resulting in a final reaction volume of 5 µL.The PCR reactions were carried out using the following cycling conditions: An initial denaturation step at 94 • C for 10 min, followed by 35 cycles of 30 s at 94 • C, 30 s at 53 • C, and 1 min at 72 • C, with a final extension step at 72 • C for 10 min.To purify the PCR products, we employed a purification protocol involving Sephadex columns (Invitrogen, Van Allen Way, North Greenbush, NY, USA).The resulting mtDNA sequences were lyophilized and subsequently resuspended in 10 µL of formamide (Hi-Di formamide AB; Applied Biosystems, Foster City, CA, USA).To obtain the sequence, both the forward-and reverse-purified PCR products were sequenced using a 3100-Avant Genetic Analyzer.

Sequences Alignment, Genetic Diversity, and Population Structure
The forward and reverse sequences of each PCR product were analyzed using SE-QUENCHER 4.2 (Gene Codes Corporation, Ann Arbor, MI, USA).Subsequently, these sequences were aligned using MUSCLE 5.0 [47].To verify the absence of stop codons, the aligned sequences were translated using MEGA 6.0.6 [48].All the obtained sequences have been deposited in the GenBank database under the accession numbers ID KT957209 to KT957275 (Table S1).
To quantify the total number and frequency of haplotypes, as well as the count of substitutions between them, we compiled a dataset consisting of all the obtained sequences, adjusting their length to 300 bp.Next, utilizing DNASP 5.10.01 [49], we calculated standard population diversity parameters, which encompassed haplotype diversity (h) and nucleotide diversity (π), at the species, subspecies, and population levels.
Additionally, we conducted Tajima's D test [50] and Fu's Fs test [51], both of which serve as measures of neutrality, using ARLEQUIN 3.11 [52].These tests allowed us to explore the possibility of past demographic expansion within the populations and subspecies.Furthermore, we employed ARLEQUIN 3.11 to analyze the geographic population structure by conducting a hierarchical analysis of molecular variance (AMOVA).The partitioning of total gene diversity was executed at three levels: Overall for the species, among the subspecies, and among the populations within each subspecies.The statistical significance of FST was assessed through 10,000 permutations.Additionally, we employed a spatial analysis of molecular variance (SAMOVA).This analysis employs simulations to identify groups of populations (k) that display geographic homogeneity while maximizing the differences between these groups.SAMOVA enables the estimation of variation between groups (FCT), between localities within each group (FSC), and between localities relative to the total sample (FST).The SAMOVA analyses were performed with 10,000 interactions for varying numbers of groups (k) ranging from 2 to 14, using the software SPADS 1.0 [53].
Furthermore, we evaluated the genetic differentiation between the subspecies and populations by estimating pairwise FST values using paired comparisons [54].This analysis incorporated haplotype frequencies and was conducted with 10,000 Markov chain steps, employed ARLEQUIN 3.11.Significance was determined following a Bonferroni correction for multiple comparisons.

Phylogenetic Reconstructions
We augmented our genetic dataset with additional Cyt-b sequences from other species within the genus Ara (A. macao, A. chloropterus, A. glaucogularis, A. ararauna, and A. severus) along with Primolius auricollis, which were accessible on GenBank.These sequences served as outgroup taxa for the phylogenetic analyses and allowed us to include fossil calibration points for the divergence time analysis [55].This augmented dataset comprised 80 terminals, including 13 outgroup taxa and 67 ingroup taxa, with a total alignment length of 582 positions.
Recent findings suggest that using multiple and distantly related outgroups can bias phylogenetic reconstructions, whereas the accuracy of phylogenetic inferences is better if the single most closely related outgroup is used [56].To assess the potential bias that the inclusion of multiple outgroups could have on the inference of phylogenetic relationships among Cyt b haplotypes of A. militaris, we also constructed a reduced dataset (with 68 terminals and 582 positions) containing all ingroup samples and a single outgroup species, A. macao, which is the closest relative to the A. militaris/A.ambiguus sister group [13,29,30].Although we recognized that the reliability of phylogenetic inferences for the A. militaris complex would increase by including its sister species, we were unable to obtain samples of A. ambiguus nor sequences of other loci that were used in previous phylogenetic studies that included both A. militaris and A. ambiguus.
To determine the optimal partitioning scheme and substitution models for both datasets (multiple-outgroups dataset and single-outgroup dataset), we initially treated each codon position as a separate partition, and we used the greedy algorithm implemented in PARTI-TIONFINDER 1.1.1[57].This analysis indicated that the multiple-outgroups dataset was optimally divided into two partitions: Partition 1 encompassed the first and second codon positions, while partition 2 comprised the third codon position within the alignment.The most appropriate substitution models identified for these partitions were HKY + G for partition 1 and GTR for partition 2. For the single-outgroup dataset, all codon positions were integrated into a single partition, evolving according to an HKY + G substitution model.
To reconstruct the phylogenetic relationships, we employed Bayesian inference (BI) and maximum likelihood (ML) approaches for both datasets.For the BI analysis, we utilized the program MrBayes 3.2 [58].For each dataset, specific substitution rates were assigned to each partition according to the PARTITIONFINDER results.Moreover, P. auricollis and A. macao were designated as outgroups for the multiple-outgroup and single-outgroup datasets, respectively.For both cases, the BI analyses comprised four independent runs, with each run extending for 10,000,000 generations and utilizing four chains.We sampled every 1000 generations to ensure a comprehensive representation of the posterior distribution.The convergence of parameters and the appropriateness of mixing in the independent runs were evaluated using Tracer 1.5.0 [59].All parameters in the analysis reached ESS values exceeding 200, signifying convergence of the chains.A definitive consensus tree was established by merging the outcomes of the four independent runs after discarding the initial 25% of the samples to ensure stability.For the ML analysis, we employed the software RAxML 7.2.8 [60].The GTRGAMMA model was separately applied to each partition for both datasets, designating the outgroups as formerly specified for the BI analyses.The reliability of nodes in the most likely tree was assessed using 1000 bootstrap replicates.Moreover, to analyze the evolutionary relationships among the identified haplotypes, we utilized the parsimony criterion implemented in Network 4.5.1.6[61] to construct a haplotype network.

Divergence Time Estimation
To estimate the divergence time between the identified lineages of the Military Macaw, we utilized the Bayesian framework within the program BEAST 2.0.1 [62].For calibration purposes, we incorporated two sources of information.Firstly, we employed a recent fossil record dating back to the late Pleistocene (approximately 1.8 million years before the present day), represented by the extinct Cuban Macaw (A. tricolor) [23].This fossil provided a minimum age limit for the genus Ara.Additionally, Tavares et al. [13,16] previously estimated the separation between the genus Ara and Primolius to have occurred 13.6 million years before the present day, with confidence intervals spanning from 9.9 to 18 million years.This calibration was established based on significant biogeographic events and contributed to the determining of divergence times among Neotropical parrots.We integrate these age constraints in the divergence time analysis.For the age of the genus Ara, we employed an offset of 1.8 Ma, which approximates the age of the A. tricolor fossil.We assigned a lognormal distribution with a mean of 3.5 and a standard deviation of 0.75.This yielded a median age estimate of 4.4 Ma, with the 97.5% Prior Credible Interval (PCI) falling below 13.3 Ma.Furthermore, we employed a normal distribution to constrain the divergence between Ara species and P. auricollis, utilizing a mean of 13.6 Ma and a standard deviation of 2.3.This distribution encompassed 97.5% of the PCI within the range of 9.09 and 18.1 Ma based on Tavares et al. [13,16].
In the divergence time analysis, we concentrated on a more compact dataset, comprising A. militaris haplotypes and a single representative from each outgroup taxa (A.macao, A. chloropterus, A. glaucogularis, A. ararauna, A. severus, and P. auricollis).This dataset was partitioned according to the previously outlined scheme.Specifically, partition 1 employed the HKY substitution model, while partition 2 used the HKY + G substitution model.We adopted a relaxed uncorrelated lognormal clock and employed a Birth Death Model as tree prior.The analysis consisted of two independent runs, each running for 100,000,000 generations, sampling parameters every 10,000 generations.We assessed the convergence and ensured satisfactory mixing of the independent runs using Tracer 1.5.0 [59].All parameters in the analysis achieved ESS values greater than 200.To generate a maximum clade credibility consensus tree, we merged the outcomes of the four independent runs and discarded 25% of the samples obtained prior to stability using TreeAnnotator 1.8.1 [58].

Contemporary and Paleo-Distribution Models: Niche Divergence Analyses
In our ecological approach, we employed ecological niche modeling (ENM) to analyze the variation in geographical and environmental conditions across the distributions of the two identified lineages of A. militaris: A. m. mexicanus lineage and A. m. militaris/A.m. bolivianus lineage.ENM offers an approximation of the Grinnellian niches occupied by each lineage, considering both contemporary and past climate conditions [63].For the execution of this analysis, we compiled occurrence data from our genetic sampling sites, along with occurrence records sourced from the Atlas of the Birds of Mexico and the Global Biodiversity Information Facility (GBIF) database [64].To account for spatial autocorrelation, we applied a filter to the occurrence records using the "spThin" R package, selecting records that were at least 5 km 2 apart [65].Consequently, the final dataset comprised spatially rarefied records, encompassing 319 occurrences for A. m. mexicanus and 299 occurrences for the South American group (A. m. militaris and A. m. bolivianus).
To define the accessibility area (referred to as "M" according to the BAM diagram) [66], we employed the unequivocal occurrence records of each lineage, which were subse-quently intersected with terrestrial ecoregions [67] and biogeographic provinces of the Neotropics [68].To characterize the ecological niche of the lineages, we utilized nine independent and uncorrelated climatic variables derived from WorldClim 1.4 [69] with a spatial resolution of approximately 5 km 2 (0.041665 • ).These climatic variables encompassed the following parameters: Mean diurnal range (Bio 02), isothermality (Bio 03), mean temperature of wettest quarter (BIO 08), annual precipitation (Bio 12), precipitation of wettest month (Bio 13), precipitation of driest month (Bio 14), precipitation seasonality (Bio 15), precipitation of warmest quarter (BIO 18), and precipitation of coldest quarter (Bio 19).The selection of these climatic variables was made based on considerations such as Pearson's correlation coefficient (r < 0.75) and Variance Inflation Factor (VIF < 10), which were evaluated utilizing the "corrplot" [70] and "usdm" [71] R libraries.
In an initial approach, we analyzed the role of ecological conditions in the differentiation among A. militaris lineages by quantifying the degree of climatic niche overlap between them [72].This enabled us to determine whether the niches of the two lineages are identical or if they are more or less similar than expected by chance.For this analysis, we employed a PCA-env approximation, a method that portrays the ecological niche in a two-dimensional space defined by the first and second principal components [73].Following this, we calculated the niche overlap between lineages, utilizing Schoener's D index.This index varies between 0 (indicating no niche overlap) and 1 (indicating complete niche overlap) as suggested by Broennimann et al. [74].To test for niche equivalency and similarity between lineages, we conducted a comparison by contrasting the empirically observed Schoener's D value with values obtained from 1000 pseudo-replicated datasets.These datasets were generated to allow for random shifts (rand type = 2, "greater" option), in accordance with Warren et al. [72] and Broennimann et al. [74].The hypothesis of niche conservatism between lineages is acceptance based on whether the observed Schoener's D value is statistically significantly higher (p < 0.05; indicating greater similarity) than the values expected from the simulated overlap.All these analyses were executed using the "ecospat v. 3.0" R library [75].
In a second approach, to understand the role of geographical drivers in the diversification of A. militaris, we generated suitability maps that represented the potential spatial-temporal geographic distribution during the Late Pleistocene climate fluctuations.To achieve this, we utilized MaxEnt 3.4.4k[76] and the "kuenm" R package [77] model for the potential environmental suitability areas for each identified lineage.First, we performed the ENM calibration protocol in the current distribution to select the best modelling parameters [78] and assessed the statistical performance of the final models [77,78].Then, we projected the determined environmental suitable conditions onto past climatic variables during the Late Pleistocene climate fluctuations, including the Mid-Holocene, Last Glacial Maximum (LGM), and Last Inter Glacial (LIG) periods.We obtained binary (presenceabsence) maps by establishing a decision threshold equivalent to the 10th percentile training presence logistic threshold.In this study, past climate data were based on three global climate circular models: CCSM4, MIROC-ESM, and MPI-ESM-P.Furthermore, we verified the reliability of our model projections by calculating the mobility-oriented parity (MOP) metric, see [79,80].
We evaluated the potential historical range extension, connectivity, and stability areas that contributed to the evolution and preservation of A. militaris lineages, e.g., [80].This was accomplished by assessing the degree of geographical matching between the potential ranges predicted for different lineages across each climatic scenario, referred to as alloprediction values.The underlying anticipation was that models would predict similar potential distribution ranges if the climatic niche of the lineages were similarly distributed and historically aligned, even in the presence of geographical barriers, e.g., [35].Our analysis led us to the identification of regions characterized by long-term climatic stability areas, often recognized as refugia.These areas were determined by overlapping the four binary maps that represent the occurrence of each lineage in all historical scenarios analyzed, e.g., [35,81].To examine whether the geographical or ecological boundaries between lineages are associated with steep environmental gradients, we utilized the linear and blob range-break tests available in the "ENMTools 1.0" R library, see [82].The scripts and input information used for all these analyses can be found at https://github.com/davidprietorres/ara_militaris_project(accesed on 25 May 2023).

Haplotype Variation and Distribution
We obtained Cyt-b sequences from 67 specimens of A. militaris.The overall count of polymorphic sites accounted for 7.04% of the entire sequence length, with parsimoniously informative sites making up just 2.8%.Among the species, we identified 41 unique Cyt-b haplotypes: 29 within A. m. mexicanus, 8 within A. m. militaris, and 4 within A. m. bolivianus (Table 2).
No shared haplotypes were found between subspecies.In particular, we identified four unique haplotypes (Hap 3, Hap 17, Hap 23) with extensive distribution within A. m. mexicanus populations.Hap 17 and Hap 23 were geographically restricted to the Sierra Madre Occidental/Sierra Madre del Sur populations.Hap 23 was present in SMOc2, SMOc3, SMOc5, and SMOc6 populations, while Hap 17 was detected in SMOc4 and SMOc5 populations.Hap 3, however, was geographically restricted to Sierra Madre Oriental populations (SMOr1 and SMOr2).For A. m. militaris populations, we identified only a single unique haplotype (Hap 40) that was widely distributed and shared between the C1 and C2 populations (Figure 1, Table 2).
Regarding the estimations of genetic diversity, in A. militaris, the haplotype diversity (h) was 0.8382, and the total nucleotide diversity (π) was 0.0063.For A. m. militaris, these values were 0.9848 for h and 0.0115 for π.In the case of A. m. bolivianus, the corresponding values were h = 0.9000 and π = 0.0097 (Table 2).As for A. m. mexicanus, the values were h = 0.7760 and π = 0.0045, indicating relatively similar levels of genetic diversity across populations of this subspecies (Table 2).Moreover, for most populations of A. m. mexicanus, no significant deviations from neutrality were observed in Fu's Fs and Tajima's D neutrality tests.The only exception was the SMOr2 and SMOc6 populations, in which negative and significant values were identified in either both tests (SMOr2: FS = −4.5029,D = −1.6238)or just one (Fs = −1.6221),implying recent expansion of these populations (Table 2).
The global AMOVA analysis revealed that most of the genetic variation (87%) was observed among populations with significant FST values (FST = 0.1304, p < 0.001).This analysis also unveiled significant differences among subspecies (FST = 0.0549, p < 0.05), suggesting that the majority of the genetic variation (66%) was observed among populations.Genetic differentiation was found between A. m. mexicanus and A. m. militaris (FST = 0.0006, p < 0.05), between A. m. mexicanus and A. m. bolivianus (FST = 0.0134, p < 0.05), and between A. m. militaris and A. m. bolivianus (FST = 0.2170, p = 0.05).In contrast, significant genetic differentiation was detected between populations (FST = 0.4612, p < 0.05) of all subspecies.For the case of A. m. mexicanus populations, the analyses indicated significant genetic differentiation between populations of the Sierra Madre Occidental/Sierra Madre del Sur and populations of the Sierra Madre Oriental (FST = 0.0803, p < 0.05).
The SAMOVA analysis identified that the A. militaris populations in this study were most effectively divided into four distinct genetic and geographic assemblages, with an FCT = 0.80, signifying 45.33% variation between the groups.The calculated FST value stood at 0.75, suggesting genetic differentiation (Table 3).According to this analysis, the A. m. mexicanus subspecies was subdivided into three assemblages: (i) Populations from the Sierra Madre Oriental (SMOr1 and SMOr2), (ii) populations from the northern of the Sierra Madre Occidental (SMOc1 and SMOc2), and (iii) populations from the central and southern portion of the Sierra Madre Occidental/Sierra Madre del Sur (SMOc3, SMOc4, SMOc5, and SMOc 6).The populations of the A. m. militaris and A. m. bolivianus subspecies constituted the fourth assemblage (Table 3).These results were reinforced by the pair-wise comparison FST analysis, which highlighted the differentiation of populations A. m. mexicanus of the Sierra Madre Oriental with those of the Sierra Madre Occidental/Sierra Madre del Sur.Additionally, differentiation was observed between the A. m. mexicanus populations and those of the South American subspecies (A. m. militaris and A. m. bolivianus) (Table 4).

Phylogenetic Reconstruction and Divergence Time Estimation
BI and ML analysis yielded comparable topologies among them and for both datasets (multiple-outgroup and single-outgroup), primarily differing in the support values (Posterior Probability [PP] and Bootstrap [BS]) assigned to the nodes (Figures 1 and S1).Therefore, in the following lines, we only describe the results of the multiple-outgroup dataset.In BI and ML analysis, the monophyly of the Ara genus and the close phylogenetic relationship between A. militaris, A. macao, and A. chloropterus were consistently recovered with relatively high support values (PP = 99% and BS = 78%).Likewise, strong support was found for the monophyly of the A. militaris species (PP = 98% and BS = 91%), as well as for clades encompassing all haplotypes of the South American subspecies, A. m. militaris and A. m. bolivianus (PP = 100% and BS = 96%).However, within this clade, only A. m. bolivianus was identified as a monophyletic lineage but with very low support (PP = 67% and BS < 50%).Furthermore, despite the unresolved phylogenetic relationships between the haplotypes of A. m. mexicanus in both analyses, we discern at least two geographically structured monophyletic groups within this North American subspecies, each with moderately high support values (PP ≥ 98% and BS ≥ 58%).The first group (G1) comprised all individuals from the Sierra Madre Oriental, while the second group (G2) exclusively consisted of individuals from the northeastern portion of the Sierra Madre Occidental (Figure 1).The rest of the haplotypes of A. m. mexicanus with unresolved phylogenetic relationships pertain to individuals from the central and southern portion of the Sierra Madre Occidental/Sierra Madre del Sur (Figure 1).
The haplotype network closely paralleled the outcomes of the BI and ML analysis, revealing a prominent phylogeographic structure in A. militaris.Within this network, we recognized two major lineages of haplotypes, distinguished by multiple mutations stapes.The first lineage encompassed all samples from the South American subspecies, wherein the haplotypes of each subspecies (A. m. militaris and A. m. bolivianus) displayed a closer relationship to one another.The second one comprised all haplotypes of the North American subspecies (A. m. mexicanus), revealing the presence of three distinct haplogroups, each uniquely confined to different geographic regions of Mexico: Sierra Madre Oriental, northern Sierra Madre Occidental, and central and southern portion of the Sierra Madre Occidental/Sierra Madre del Sur (Figure 1).Notably, we also observed that most of the haplotypes of A. m. mexicanus slightly differed from the most common and widely distributed haplotypes (H23, H29, H17, and H13) in this subspecies, whereas a larger number of differences was observed among the haplotypes of the South American subspecies.
Our divergence time analysis suggests that the separation between A. militaris and A. macao, as well as A. chloropterus, likely occurred approximately 2.35 Ma BP (95% Confidence Interval, CI, ranging from 4.87 to 0.51 Ma BP).Furthermore, this analysis also indicates that the mitochondrial lineages divergence within A. militaris was likely initiated around 1.56 Ma, during the Pleistocene (Figure 1).However, it is important to note that the 95% confidence interval for this estimated mean divergence time was relatively large ranging from 3.89 to 0.26 Ma before the present day.

Niche Divergence Analyses and Paleo-Distributions Models
Our models showed a good fit for the current distribution for both the Mexican (partial-ROC test = 1.50; omission rate = 13.70%;AICc = 4,764.19)and the South American lineages (partial-ROC test = 1.46; omission rate = 14.6%;AICc = 5098.47).In both instances, the results indicated strong model performance, indicating that the models surpassed random outcomes and were therefore statistically descriptive of the Grinnellian niche of each lineage.The Jackknife test and variable contribution values within MaxEnt highlighted the most influential variables in the model construction for the A. m. mexicanus lineage, which included BIO 12 (34.8%),Bio 15 (16.0%),BIO 08 (15.6%), and Bio 03 (13.3%).In the case of the South American lineage (A. m. militaris and A. m. bolivianus), the most important variables were BIO 08 (38.9%),Bio 13 (23.3%),and BIO 19 (10.4%).
In contrast, the occurrence of density surfaces (dark spots in Figure 2) in environmental space, as determined by PCAenv, distinctly portrayed the differentiation in positioning between the South and North American lineages.Results from the PCA-env revealed that 72.50% of the environmental variation is explained by the two first principal components (PC1 = 59.75% and PC2 = 12.75%).Observed Schoeners' D values were 0.14 and statistically non-significant in both niche equivalency (p = 0.99) and similarity tests (A. m. mexicanus-South American group: p = 0.12 vs.South American group-A.m. mexicanus: p = 0.21).These results suggest that the observed niche overlap was no more similar/equivalent than expected by chance (Figure 2), thus rejecting the hypothesis of niche conservatism between lineages (that is, lineages occupy different niches).Finally, both the linear and blob range-break tests supported the idea that significant environmental differences (p < 0.05) separate the lineages.0.21).These results suggest that the observed niche overlap was no more similar/e lent than expected by chance (Figure 2), thus rejecting the hypothesis of niche con tism between lineages (that is, lineages occupy different niches).Finally, both the and blob range-break tests supported the idea that significant environmental diffe (p < 0.05) separate the lineages.We only present values for the D metrics for both tests.In all plots, the arrow with a red diamond represents the observed similarity value between niches; the gray columns represent 1000 randomly simulated expected values.The p-value is shown in each plot (n.s.= non-significant).
Furthermore, results from the paleodistribution models (Figure 3) demonstrated that the suitable areas' distribution for the two lineages was markedly geographically differentiated.Notably, the allopredictions (i.e., predicted distributions of the Mexican lineage into the range of the South American lineage and vice versa) exhibited lower overlap values across all four scenarios (5.6% for the Recent, 6.7% for the mid-Holocene, 2.5% for the LGM, and 9.0% for the LIG).This suggests that environmental differentiation has been upheld since at least the Late Pleistocene.Overall, we observed that the South American lineage demonstrated an increase (~66%) in highly suitable areas during the LGM in comparison to the present, whereas the Mexican lineage displayed relatively smaller suitable areas (approximately 57%) that were also more fragmented.Moreover, the projected climatic stability areas for each lineage (crossed areas within the current distribution map, Figure 3) covered 38.62% of the presently recognized distribution for the Mexican lineage and 52.56% for the South American lineage.However, it is important to highlight that models predicted a higher surface of climatic stability areas for South American lineage across Central America and northwestern Mexico (Figure 3).
areas (approximately 57%) that were also more fragmented.Moreover, the projected climatic stability areas for each lineage (crossed areas within the current distribution map, Figure 3) covered 38.62% of the presently recognized distribution for the Mexican lineage and 52.56% for the South American lineage.However, it is important to highlight that models predicted a higher surface of climatic stability areas for South American lineage across Central America and northwestern Mexico (Figure 3).

Genetic Diversity and Structure
We detected relatively high levels of haplotype diversity but low nucleotide diversity in the Military Macaw, which is consistent with the high frequency of private haplotypes and limited genetic differentiation [83].This genetic diversity profile is akin to that observed in the endangered species Hyacinth Macaw (Anodorhynchus hyacinthinus) [84] and in psittacidae species classified of least concern such as the Scarlet Macaw (A. macao) [85] and Burrowing Parakeet (Cyanoliseus patagonus) [17].
Nonetheless, the overall genetic diversity we observed in the Military Macaw surpasses that of other threatened bird species [86][87][88].Despite experiencing population declines and increased fragmentation due to habitat loss [24,25,89], the relatively elevated levels of genetic diversity unveiled in this study suggest historically large population sizes for the species.Comparable scenarios have been proposed for other endangered bird species, including the Andean Condor (Vultur gryphus) [90], Whooping Crane (Grus americana) [43], Greater Prairie-Chicken (Tympanuchus cupido) [91,92], the Spotted Owl (Strix occidentalis) [93], and the Resplendent Quetzal (Pharomachrus mocinno) [94].Furthermore, we observed higher levels of genetic diversity in A. m. militaris and A. m. bolivianus in comparison to A. m. mexicanus, implying that the South American populations historically maintained larger population sizes than the North American subspecies of the Military Macaw.
Our findings highlight the existence of two major populations within the Military Macaw.The Mexican population (subspecies A. m. mexicanus) shows significant genetic differentiation from the South American populations (A. m. militaris + A. m. bolivianus), with no shared haplotypes between them.This partitioning can be attributed to long-term population isolation, as observed in other macaw species such as Ara ararauna [8].In contrast, for the populations of the subspecies A. m. bolivianus and A. m. militaris, no significant FST values were obtained.This implies that these taxa continue to experience gene flow, which is plausible given their geographical proximity.Nevertheless, this pattern can also be explained by incomplete lineage sorting within the mitochondrial genome of these subspecies [32].
We observed genetic differentiation among populations inhabiting the Sierra Madre Occidental/Sierra Madre del Sur and Sierra Madre Oriental in A. m. mexicanus.A previous study conducted in the same Mexican localities using nuclear markers (microsatellites) [37] reported limited gene flow between these regions, likely attributed to geographic barriers such as mountain ranges (Trans-Mexican Volcanic Belt, Sierra Madre Occidental, and the Sierra Madre Oriental, as well as the Central Mexican Plateau) [95,96].Additionally, our study revealed further sub-structuring among populations in the northwestern Sierra Madre Occidental/Sierra Madre del Sur region.This pattern was not evident in previous microsatellite data [37] but can be explained by sex-biased dispersal [6,8,96,97].For instance, if dispersal is predominantly male-biased while females tend to exhibit philopatry, nuclear markers are more likely to detect higher levels of gene flow among populations compared to mitochondrial markers [8].Furthermore, in such a scenario, a significant proportion of localized mitochondrial haplotypes can be expected [8,98,99], mirroring the observed in our present study.

Phylogeography, Niche Divergence, and Paleo-Distributions Models
Our phylogenetic analysis revealed significant differentiation between the haplotypes of A. m. mexicanus and two South American subspecies (A. m. militaris and A. m. boliviaunus), providing support for the existence of two distinct lineages within A. militaris.This outcome aligns with previous phylogenetic inferences based on different mitochondrial genes [9] or on the whole mitochondrial genome of the species [32].It is noteworthy that the populations of A. militaris from North America and South America are entirely allopatric.Furthermore, our results suggest the historically limited movement of mitochondrial lineages between these regions, as indicated by the absence of shared haplotypes among populations [8].Additionally, our estimates of divergence time estimations indicate that the genetic divergence between these two populations likely occurred during the Pleistocene.
The repeated instances of discontinuous distributions observed in A. militaris can potentially be attributed to two distinct processes: Long-distance dispersal or local extinction of a once widely distributed ancestor [100].However, it is important to note that long-distance dispersal is considered less likely in non-migratory avian taxa [100].Furthermore, although individuals of A. militaris can fly distances of around 24 km per day, these movements are seasonal, altitudinal, and nomadic, primarily in response to fluctuations in food availability [37,101,102].Additionally, there are no records of any species within the genus Ara undertaking latitudinal migratory movements [25,37].These various lines of evidence collectively suggest that long-distance migration is an unlikely explanation for the discontinuous distribution observed in A. militaris.
On the other hand, extinction has been widely considered one of the most plausible explanations for the broad distributional gaps observed in non-migratory avian taxa, e.g., [100,103].The disjunct distribution pattern of A. militaris could potentially be attributed to local extinctions that occurred relatively recently, possibly associated with the effects of Pleistocene climatic changes around 1.6 million years ago.In Central America, the Pleistocene climatic cycling between glacial and interglacial periods led to significant changes in the distribution and population sizes of various plant and animal species, as evidenced by fossil and genetic data spanning the Quaternary period (last 2.6 million years).These data demonstrate fluctuations in the composition and distribution of species currently found in dry forests [104,105].During glacial periods, cooler and drier conditions prevailed in the Neotropics, facilitating the expansion and connectivity of dry forests across the Americas [106,107], including suitable habitats for species such as A. militaris.
In contrast, during interglacial periods, characterized by the conditions experienced today, warm and humid environments prevail.This can potentially lead to the disruption of the dry forest corridor due to the expansion of humid forests in the Neotropics [108,109].Such expansion could significantly reduce the extent of dry forests, ultimately resulting in the local extinction of A. militaris populations in Central America and causing a separation between populations in North America and South America.The results of our paleodistribution models align with this extinction scenario in two main ways.First, for the South American lineage, a larger distribution was predicted during the last glacial maximum compared to the three different interglacial periods assessed (current, Holocene, and LIG).Second, although historically stable areas are predicted in Central America for this lineage, they are relatively small and widely fragmented, potentially harboring small populations that are more susceptible to local extinction events.
In addition to historical habitat reductions, the extinction of A. militaris in Central America could have been influenced by negative interactions with its sister species, A. ambiguus.This species primarily inhabits Central America and possesses a larger body size compared to the Military Macaw.While A. ambiguus is commonly associated with humid and lowland forests [110], these environmental conditions could provide a potential opportunity for A. ambiguus to outcompete A. militaris, potentially contributing to the latter's extinction in Central America.However, this hypothesis warrants further investigation and should be explored more comprehensively.
In Mexico, we have identified two geographically structured monophyletic groups of A. m. mexicanus.One group is restricted to northeast Mexico in the Sierra Madre Oriental, while the other is confined to the northwestern portion of the Sierra Madre Occidental/Sierra Madre del Sur.These findings align with previous molecular studies on the species using mitochondrial [31] and nuclear [37] markers.However, our results also suggest additional genetic structure among the populations in western Mexico.The phylogeographic and genetic structure revealed in this study indicates historical isolation among the populations identified in Mexico.This isolation could be attributed to the presence of geographic barriers [37,111] as well as the responses of A. m. mexicanus populations to Quaternary climatic fluctuations [112].Indeed, our paleodistribution models for the North American lineage of A. militaris indicate historically small and fragmented stable areas in Mexico.Additionally, we found that the distribution of this lineage has undergone changes during the climatic dynamics of the last glacial and interglacial periods, although in a different manner than observed for the South American lineage.Specifically, our models predict a smaller and more fragmented distribution during the glacial periods compared to the interglacial times.It is worth noting that palynological studies have indicated that climatic conditions in Mexico during interglacial or glacial times were not uniformly dry or wet, respectively [113][114][115].Thus, during glacial periods, the dry forest biome could contract in response to more humid conditions, whereas it could expand during interglacial periods due to drier conditions.Similar responses to Quaternary climatic conditions, as observed in A. m. mexicanus, have also been identified in other dry forest species.For example, similar patterns have been observed in the Long-billed hermit (Phaethornis longirostris) [116], Red-eyed Vireo (Vireo hypochryseus) [33], Russet-crowned Motmot (Momotus mexicanus), Golden-cheeked Woodpecker (Melanerpes chrysogenys), Orange-breasted Bunting (Passerina leclancherii) [34], and Squirrel Cuckoo (Piaya cayana) [35].
If populations of A. m. mexicanus contracted during glacial periods and subsequently expanded from refugia during interglacial periods, we would expect to find evidence of a recent population expansion in the species.In line with this expectation, neutrality tests indicate a recent demographic expansion in at least two Mexican populations.Additionally, the star-like relationships observed in the haplotype network analysis for the Mexican populations are commonly observed in other taxa that have undergone recent demographic expansions, including invasive species [117,118].

Conservation and Taxonomic Implications
The destruction and fragmentation of natural forests due to increased human activities have had a significant impact on the natural habitat of A. militaris.Furthermore, the capture and illegal trade of wild macaw individuals have played a major role in the decline of wild populations [21,22,37].In this context, conservation genetics has played a crucial role in developing appropriate conservation and management strategies [119].One approach is the identification of conservation units, such as Evolutionarily Significant Units (ESUs) and Management Units (MUs) [119][120][121].ESUs and MUs help to identify patterns in genetic variation that result from evolutionary history, enabling us to make informed decisions for the conservation of genetic diversity within a species [121].In fact, the phylogeographic patterns of the Military Macaw identified in this study could contribute to defining conservation units, especially considering that most of the individuals included in the study were from breeding populations.
According to the Moritz model [36], the populations of Military Macaw in South America (A.m. militaris and A. m. bolivianus) should be considered a single ESU.This is further supported by the results of the ecological niche divergence analysis, which indicate that the environmental niches occupied by the North American and South American lineages are distinct.Therefore, the conservation efforts for these lineages should be separated.Within the ESU in South America, two subgroups can be distinguished corresponding to the subspecies A. m. militaris and A. m. bolivianus.These subgroups can be recognized as independent MUs because they possess unique and non-shared haplotypes [119,120,122].
For the Mexican populations of the Military Macaw, it is appropriate to recognize at least two ESUs.The first ESU corresponds to the populations in the Sierra Madre Oriental, as their genetic differentiation has been consistently observed using both mitochondrial markers [9,32] and nuclear [37] markers.The second ESU corresponds to the populations in the Sierra Madre Occidental/Sierra Madre del Sur.Within this ESU, two MUs can be considered: (1) The northern portion of the Sierra Madre Occidental and (2) the central and southern portion of the Sierra Madre Occidental/Sierra Madre del Sur.These regions show distinct genetic structures and differentiation.However, it is worth noting that Rivera-Ortiz et al. [37] considered the northwestern populations as a single MU.
Overall, our findings are largely consistent with previous phylogenetic studies on the Military Macaw, which were based on a larger number of mitochondrial loci [31], or the entire mitochondrial genome [32].These studies identified two major lineages within the species: The North American subspecies (A. m. mexicanus) and the South American subspecies (A. m. militaris and A. m. bolivianus) [31,32].These results suggest a historical isolation between the ancestral populations of the military macaw in North and South America, a separation that still persists in its discontinuous geographic distribution.Additionally, we have identified significant ecological niche differentiation between these two lineages.Taken together, these findings suggest that A. m. mexicanus can be recognized as a distinct species (Ara mexicanus) separated from the Military Macaw in South America, which can improve the conservation priorities for both taxa [94,120].Nevertheless, genetic information from nuclear regions, as well as morphological measurements and ecological aspects, are indispensable to further support the taxonomic status of these two lineages [32,35].
Finally, we recognize that our inferences about the phylogeographic history of A. militaris can be limited by the use of short sequences of a single locus and the omission of its sister species (A.ambiguus).Although, our results were largely consistent among them, as well as with previous studies including a larger number of mitochondrial loci and samples of A. ambiguus [29,30], we also consider that further studies involving a broader set of independent markers and populations (of both A. militaris in South America and

Figure 1 .
Figure 1.Geographic distribution and phylogenetics of the Military Macaw.(A) mtDNA phylogenetic tree based on a Cyt-b gene (>300 bp) concatenated dataset.Numbers in branches indicate nodal support based on Bayesian posterior probabilities and strict consensus of maximum-likelihood.Numbers above branches indicate divergence time in red.(B) Geographic distribution in the natural range and distribution of Cyt-b haplotypes in the subspecies (A. m. militaris, A. m. bolivianus, and A. m. mexicanus).To make reading the figure easier, the haplotypes have been grouped according to position in median-joining network of Figure 1C.The black dots represent the sample populations in the distribution range of Military Macaw.Dotted lines depict both the accessible (or M) and projection areas for the models and the niche divergence analyses for the two evolutionary lineages herein assessed.Brown shading indicates areas at least 1000 m above sea level (m.a.s.l).(C) Haplotype network.Colors correspond to the vertical bar in the phylogenetic tree.The size of each circle is proportional to the corresponding haplotype frequency.Red circles inside the branches indicate missing intermediates; the size of each branch is proportional to the mutational steps.

Figure 1 .
Figure 1.Geographic distribution and phylogenetics of the Military Macaw.(A) mtDNA phylogenetic tree based on a Cyt-b gene (>300 bp) concatenated dataset.Numbers in branches indicate nodal support based on Bayesian posterior probabilities and strict consensus of maximum-likelihood.Numbers above branches indicate divergence time in red.(B) Geographic distribution in the natural range and distribution of Cyt-b haplotypes in the subspecies (A. m. militaris, A. m. bolivianus, and A. m. mexicanus).To make reading the figure easier, the haplotypes have been grouped according to position in median-joining network of Figure 1C.The black dots represent the sample populations in the distribution range of Military Macaw.Dotted lines depict both the accessible (or M) and projection areas for the models and the niche divergence analyses for the two evolutionary lineages herein assessed.Brown shading indicates areas at least 1000 m above sea level (m.a.s.l).(C) Haplotype network.Colors correspond to the vertical bar in the phylogenetic tree.The size of each circle is proportional to the corresponding haplotype frequency.Red circles inside the branches indicate missing intermediates; the size of each branch is proportional to the mutational steps.

Figure 2 .
Figure 2. Equivalence and similarity test in the environmental space for the two herein-identified lineages of Military Macaw.(A) A Principal Component Analysis (PCA) of ecological niche for the lineages and the corresponding values of variables contribution.(B) Dark shading in the upper right plots depicts the density of the occurrences of each lineage by cell in the ecological space; green shading is the density of the occurrences of A. m. militaris, while the red color is the density of the occurrences of A. m. mexicanus.The solid and dashed lines indicate 100 and 50% of the available (background) environment, respectively.(C) Plot of the equivalence test comparing the ecological niche of two lineages.(D) Plots of the similarity test comparing the ecological niche of the two lineages in both directions (A. m. mexicanus vs. A. m. militaris/A.m. bolivianus and vice versa).We only present values for the D metrics for both tests.In all plots, the arrow with a red diamond represents the observed similarity value between niches; the gray columns represent 1000 randomly simulated expected values.The p-value is shown in each plot (n.s.= non-significant).

Figure 3 .
Figure 3. Ecological niche model projected onto the geographic areas for the two herein identified lineages of Military Macaw: A. m. mexicanus (blue) and A. m. militaris/A.m. bolivianus (orange).Upper maps depict suitability areas for each lineage under recent climatic conditions; crossed lines and dotted red circles depict regions of historical climatic stability.Lower maps depict the predicted regions of historical distribution for climatic niches across four time periods, and the aloprediction areas (darker blue and darker orange) between lineages across the Late Pleistocene (LIG and LGM), the Holocene and the Recent.

Table 1 .
Sample sites of the Military Macaw (Ara militaris) in North America and South America.

(Subspecies) Locality (State/Department) Abbreviation
+ Samples obtained from museum collections.º Samples obtained directly in the field.

Table 2 .
Mitochondrial DNA diversity indices for A. militaris populations, including nucleotide (π) and haplotype (h) diversities, as well as Tajima s D, and Fu s Fs neutrality test indexes.Neutrality statistics Tajima's D, and Fu's Fs are bold when significant at a = 0.05.Values ± SD, * Extirpated population.

Table 4 .
Pair-wise comparison of FST values between the 12 populations of A. militaris.