Falagonia mexicana is an aleocharine distributed from northern Mexico to Guatemala and El Salvador. It is associated with Atta mexicana ants and lives within their piles of waste or external debris. The phylogeography and historical demography of 18 populations from Mexico, Guatemala, and El Salvador were studied. The data set encompasses a 472 bp fragment of the COI. Results suggest that F. mexicana was originated during Middle Pliocene (ca. 0.5 Mya), starting its diversification at the Upper Pleistocene and Holocene. Populations were recovered forming at least four main lineages, with a significant phylogeographic structure. Evidence of contemporary restricted gene flow was found among populations. The historical demography suggests that the geographic structure is due to recent physical barriers (e.g., Isthmus of Tehuantepec) rather than ancient geological events. Also, recent geological and volcanic events in the east of the Trans-Mexican Volcanic Belt and the Sierra Madre Oriental might be responsible for the restricted gene flow among populations. Skyline-plot analyses suggested that a demographic expansion event took place at the end of the Late Quaternary glacial-interglacial cycles.

Cytochrome oxidase I, “false Lomechusini”, Mesoamerica, mtDNA

Falagonia mexicana Sharp, 1883 is a species of Coleoptera that belongs to the subfamily Aleocharinae (Staphylinidae) (Bouchard et al. 2011). Particularly, the species is placed within Lomechusini, a tribe distinctive for having a 4–5–5 tarsal formula, metasternal process typically longer than mesosternal process, maxillar galea and lacinia moderately to considerably elongated, mesocoxae moderately to very widely separated by broad meso- and metaventral processes (Seevers 1965, 1978), among other characters. Several Lomechusini species are associated with termites and ants (Navarrete-Heredia et al. 2002), as is the case of F. mexicana. However, some recent studies have suggested that Lomechusini is not a monophyletic group and two clades have been proposed: the “true Lomechusini” whose species show a Holarctic distribution, and the “false Lomechusini” distributed in tropical America (Elven et al. 2010, 2012). Based on molecular data F. mexicana has been recovered as part of the “false Lomechusini” (unpublished data).

The genus Falagonia Sharp, 1883 included two species: F. mexicana and F. crassiventris Sharp, 1883. Currently, F. crassiventris is located in the genus Pseudofalagonia Santiago-Jiménez, 2010; therefore, Falagonia is now a monotypic genus (Santiago-Jiménez 2010; Santiago-Jiménez and Espinosa de los Monteros 2016). At a phylogenetic analysis, Falagonia was supported by seven autapomorphies: 1) epipharynx with a spinose process on the apico-medial margin; 2) with a pair of setae on the sensory area of ligula; 3) with simple punctures on elytral disc; 4) without tuberculate punctures on elytral disc; 5) tergite IV of male with a protuberance on each side of midline; 6) tergite VIII of male without lobes on posterior margin; and 7) tergite VIII of female without lobes on posterior margin. It is important to highlight that the specimens of F. mexicana are very homogeneous, except for some variation in size (4–7 mm); in color (from reddish to pale yellow in the tenerals); in the length of the α-sensilla of the epipharynx, as well as the length of the spines of the spinose process on the apico-medial margin of the epipharynx. Moreover, they present secondary sexual dimorphism. Males have a protuberance on each side of the midline of tergite IV, and a protuberance on the midline of tergite VIII.

From Mexico to El Salvador, F. mexicana has been recorded in a wide altitude range (from 50 to 2,000 m a.s.l.). This species has been recorded in diverse types of vegetation (Santiago-Jiménez 2010). According to Challenger and Soberón (2008) those habitats correspond to tropical evergreen forests, tropical deciduous forests, temperate coniferous and broadleaf forests, cloud forests, and scrub xerophytes. These types of vegetation in Mexico tend to present warm humid climates (Af, Am, Aw, BSh) to temperate climates (Cfb, Cwa, Cwb). Records in pine and oak forests do not exceed 2,000 m a.s.l. Thus, in general, they prefer warm humid climates with dry winter (Rzedowski 1994).

In Mexico, F. mexicana is found in the states of Colima, Guanajuato, Guerrero, Hidalgo, Jalisco, Michoacan, Morelos, Oaxaca, Queretaro, San Luis Potosi, Sinaloa, Sonora, Tamaulipas, and Veracruz (Santiago-Jiménez 2010). This distribution corresponds to the biogeographic provinces of the Pacific, Trans-Mexican Volcanic Belt, Sierra Madre Oriental, Sierra Madre del Sur, Soconusco, and Tierras Altas de Chiapas. These regions are associated with physiographic provinces that show heterogeneous historical processes and dates of origin (e.g., de Cserna 1989; Mastretta-Yanes et al. 2015 for summaries). For instance, the most recent formation corresponds to the Trans-Mexican Volcanic Belt from the Miocene to the present day, with some stratovolcanoes forming in during the Pleistocene (Ferrari et al. 2012).

Falagonia mexicana lives in the debris produced by Atta mexicana (Smith, 1858) ants. The reddish color F. mexicana allows it to blend in with the debris. Apparently, despite having developed wings, this species does not fly. Flight intercept traps have been placed near the debris and we have not been able to catch them. Instead, by placing pit-fall traps in the vicinity of the detritus we have obtained some specimens. The nests of A. mexicana with its associated debris seems to be islands in the geographic distribution of F. mexicana, due to their presumed inability to fly and limited long-distance movement. Therefore, geographical barriers have a great effect in reducing gene flow between its populations, as well as restrictions on its geographical distribution. Physical barriers can lead to genetic isolation when they remain for an extended period of time, which is essential for lineage divergence (Saunders et al. 1991; Avise 2000).

Apparently, both A. mexicana and F. mexicana are incapable of surviving in areas above 2,300 m a.s.l. In the case of F. mexicana, its septentrional distribution is limited by the Tropic of Cancer, as is the case with other insect groups (e.g., Passalidae, Reyes-Castillo 2000). However, A. mexicana does reach the southern United States. It has been seen in the field that F. mexicana prefers debris that is under vegetal cover, with more shade and humidity. Usually, F. mexicana is found in A. mexicana’s detritus, but that does not mean that every A. mexicana nest has F. mexicana associated with it, especially when the detritus is exposed to desiccation.

Recently, some studies have been carried out to understand the biogeographical patterns of the biota in Middle America (e.g., Marshall and Liebherr 2000; Daza et al. 2010; Ornelas et al. 2013). With an average altitude of 200 m a.s.l. the Isthmus of Tehuantepec has been identified as an effective barrier for gene flow (Barrier et al. 1998). Likewise, the Trans-Mexican Volcanic Belt represents a biogeographical barrier (Marshall and Liebherr 2000) and a center for lineage diversification (Bryson et al. 2018). The use of coalescent and geographic structure models based on molecular data can be useful to test historical demographic scenarios where hypothesis of vicariance and dispersion can be contrasted (Knowles and Carstens 2007; Daza et al. 2010). Therefore, F. mexicana is an interesting model to study the relationship between genealogy and geography. The objectives of the present study are to estimate the date of the appearance of Falagonia mexicana, to evaluate the phylogenetic structure of F. mexicana, and to correlate the historic events that explain the divergence processes within the lineage.

The sampling was carried out in 18 localities from El Salvador, Guatemala, and Mexico (Fig. 1, Table 1). The collection sites were selected based on previous records and potential suitable sites. Potential localities were identified based on their environmental and geographic conditions. This sampling covers most of the species’ distribution. Field collections were carried out by direct collection on debris of A. mexicana using an aspirator. After its collection, the specimens were placed in 96% ethanol and stored at -20 °C. The material was deposited at the Phylogenetic Systematics Laboratory (INECOL, Xalapa).

Figure 1. 

Distribution map showing the localities where individuals of Falagonia mexicana were sampled. Pie charts represent haplotypes sampled in each locality. Color on map indicates altitudinal changes every 500 m a.s.l. The polygons represent the main biogeographical regions: Sierra Madre del Sur (black), Sierra Madre Occidental (purple), Sierra Madre Oriental (green), Tierras Altas de Chiapas (blue), and Trans-Mexican Volcanic Belt (red). The white line represents the Isthmus of Tehuantepec.

From total genomic DNA, we obtained the sequence for a fragment of 472 base pairs (bp) of the Cytochrome Oxidase subunit I (COI) mitochondrial gene. This fragment is located between positions 250 and 720 of the COI gene. This gene was selected since it presents an adequate mutation rate for this kind of studies (Avise 2000). Sequences present extensive intraspecific polymorphisms, providing information between individuals. In addition, this molecular marker has been a useful tool for species’ identification (Avise 2000). DNA extraction was performed using the kit Tissue & Insect DNA MicroPrep (Zymo Research Corp., Irvine, CA), following the protocol suggested by the manufacturer. DNA was extracted for amplification from 139 specimens. Amplification was performed using the S1718 (Normark 1996) and NANCY (Simon et al. 1994) primers designed for Coleoptera. PCR mix consisted of 5 µL of DNA, 3 µL of 5X Buffer, 3 µL of MgCl2 (25mM), 3 µL of dNTP’s (8mM), 3 µL of each primer (10mM), 0.2 µL of Taq polymerase (5 U / µL) and 9.8 µL of ddH2O adding to a final volume of 30 µL. The PCR protocol consisted of an initial denaturation for 180 s at 94 °C; 35 cycles of 30 s at 94 °C as denaturing temperature, 30 s at 45 °C for primer annealing and 90 s at 74 °C for extension; followed by one final extension cycle of 300 s at 72 °C. PCR products were visualized in a 1% agarose gel stained with ethidium bromide and purified using the GeneJET Genomic DNA Purification Kit (Thermo Fisher Scientific Inc., Waltham, MA). The PCR samples were sent to the Macrogen company (South Korea) for sequencing. Sequences were assembled with the software SEQUENCHER v. 5.2.4 (Gene Codes Corporation). The alignment of the sequences was trivial; therefore, it was done manually using the software MEGA v. 7 (Kumar et al. 2016). All the sequences acquired for the first time were deposited in GenBank (OQ608846–OQ608984).

The identity of the sequences was corroborated by translating them into amino acids, verifying the absence of nonsense codons or intermediate stops, and by a BLAST search. Using the software PAUP v. 4.0 (Swofford 2002), the nucleotide proportion of each sequence was estimated, as well as their standard deviation and nucleotide bias. In turn, the first, second, and third positions were identified, and the same descriptors were calculated. The proportion of transitions and transversions were also calculated.

To identify genetically congruent geographic regions, a spatial analysis of molecular variance (SAMOVA) was conducted. This was performed with the aid of the software SAMOVA v. 1 (Dupanloup et al. 2002), using 100 banding steps to obtain the fixation index F_CT (Excoffier et al. 1992). Likewise, a non-hierarchical analysis of molecular variance (AMOVA) was carried out including every population, followed by an analysis between the statistically significant regions identified by the SAMOVA. This was conducted with 1000 permutations using the software ARLEQUIN v. 3.5 (Excoffier and Lischer 2010). A Mantel test was performed, based on the combined data set, to detect possible isolation by distance. Genetic variation within and between populations was evaluated using descriptive empirical values, such as number of haplotypes (H), haplotype diversity (Hd), nucleotide diversity (π) (Nei 1987), and segregated sites (S) (Nei and Kumar 2000). These values were estimated using the software DNASP v. 5.10 (Librado and Rozas 2009). This analysis was done for each population, as well as for each region inferred by the SAMOVA. The haplotype network was made with the software NETWORK v. 4.5 (Bandelt et al. 1999) using the Median Joining algorithm with a weight of 10 for all characters and an epsilon of 10 as suggested by the authors. Similarly, the Maximum Parsimony post-processing tool was used to eliminate all the superfluous intermediate haplotypes and the links that were not in the shortest network.

The genealogy of the sequences was estimated with a Bayesian inference analysis, using the software MRBAYES v. 3.2.3 (Ronquist et al. 2012). The best-fit nucleotide substitution model was selected using the Akaike information criterion through the software JMODELTEST v. 2.02 (Posada 2008). A homologous sequence of Pseudofalagonia sp. (acc. num. OL764948.1), a species closely related to F. mexicana, was used as outgroup. The Markov Monte Carlo Chains (MCMC) were run for 10 million generations sampling trees and parameters every 1000 generations. The consensus tree was obtained with their respective posterior probabilities after discarding the initial 25% of the stored trees (Ronquist et al. 2012). Likewise, an inference based on the Maximum Likelihood (ML) optimality criterion was carried out using the software RAXML v. 8 (Stamatakis 2014). The robustness of the nodes recovered in the ML tree was evaluated with 500 bootstrap replications.

To date and infer possible isolation events in F. mexicana populations, we used the software BEAST v. 1.10.4 (Drummond et al. 2012). To optimize tree search, the topology inferred from the Maximum Likelihood analysis was used as the tree prior. The ML tree was used instead of the Bayesian one because BEAST only accepts fully resolved trees as priors. To estimate the divergence times, we used a relaxed uncorrelated clock model.

A nucleotide substitution rate of 6.311 (10^-2 subs/site/My/lineage) for the COI was used to date the genealogy of Falagonia. This rate between taxa is based on the substitution rate estimated by Pons et al. (2010) for the COI gene of the coleopteran mitochondrial genome. Alternatively, a more conserve substitution rate of 2.0 (10^-2 subs/site/My/lineage) was used (chronogram available as Suppl. material 1). This rate closely matches the standard mitochondrial arthropod clock reported by Brower (1994), and the rates calculated in studies that compared closely related species of Coleoptera (e.g., Leys et al. 2003; Balke et al. 2009; Ribera et al. 2010). The BEAST MCMC ran for 10 million generations, sampling trees and parameters every 1000 generations. The software TRACER v. 1.6 (Rambaut et al. 2014) was used to evaluate stationarity of the MCMC parameters, sample sizes (ESSs > 200) and posterior intervals spanning the 95% of the highest posterior density. This analysis was repeated four times and the trees resulting from each run were combined using the software LOGCOMBINER v. 1.10.4 (Drummond et al. 2012), after applying an initial burn-in of 25%. The evaluated nodes were those that scored a posterior probability value higher than 0.9. The inferred chronogram was viewed and edited with the software FIGTREE v. 1.4.2 (Rambaut 2008). To determine the historical demography of the populations, we conducted a neutrality test with Fu’s F statistics (Fu 1997). This was compared with the generalized skyline-plot analysis inferred with BEAST. The input parameters, as well as the molecular clock model, were identical to the ones used for the coalescence analysis for node dating.

The alignment of the 139 COI sequences showed 417 constant positions (88%). The 55 remaining characters correspond to variable sites (See Suppl. material 2), of which 39 of them are potentially phylogenetically informative, while the other 16 are singletons. Sequences present a 0.27 nucleotide bias and the following nucleotide percentages: A 29.8%, C 17.2%, G 14.9%, and T 38.1%. Regarding codon positions, first positions show less nucleotide bias (0.11) and the following nucleotide proportion: A 25%, C 17.6%, G 25.0%, and T 31.2%. Second positions present an intermediate bias (0.28) and the following nucleotide percentages: A 13.3%, C 28.8%, G 17.2%, and T 40.7%. Lastly, third positions show the highest nucleotide bias (0.64) and this nucleotide composition: A 50.4%, C 17.2%, G 14.9%, and T 38.1%. Regarding base substitution mutations we found a total of 48 transitions and 30 transversions, resulting in a Ts:Tv = 1.6:1 ratio. Mutation transitions were as follows: 27 A-G and 21 C-T; while transversions: 10 A-T, 10 T-G, 7 A-C and 3 C-G.

The haplotype network (Fig. 2) contains 60 haplotypes distributed among the 18 localities. Most of the localities exhibited more than one haplotype. Two main haplotypes (i.e., H1, H46) connected by eleven mutations were recovered. Most of the Sierra Madre Oriental populations had haplotype H1 individuals; while haplotype H46 is shared by individuals from the Sierra Madre del Sur, and one population from the Valles Centrales de Oaxaca. Haplotypes H6 and H20 are shared between the Sierra Madre Oriental and Tierras Altas de Chiapas regions. Haplotypes H2 and H24 collected in Neblinas (Queretaro state) were recovered in the network as the most closely related to the haplotypes from Sierra Madre Oriental and Tierras Altas de Chiapas, respectively. Haplotype H19 from Lachiguiri (Oaxaca state) is an intermediate haplotype among those of Sierra Madre Oriental, Tierras Altas de Chiapas, and Sierra Madre del Sur in Valles Centrales de Oaxaca. Individuals from Tlalixtac and Yautepec (Valles Centrales de Oaxaca) shared Haplotype H39. Haplotypes H41, H42, and H43 are endemic to Flor de Chiapas. Haplotypes H38 and H40 are endemic to Tlalixtac. Haplotypes H45, H47, H48 and H49 are endemic to Candelaria Loxicha (Oaxaca state), and haplotypes H50 to H53 to Yautepec. Haplotypes H54 to H60 were found in individuals from the Sierra Madre del Sur (Acahuizotla, Guerrero state) and the Trans-Mexican Volcanic Belt (Tepexco, Puebla state) regions. Global and local haplotype frequencies can be consulted in Suppl. material 3.

Figure 2. 

Haplotype network of Falagonia mexicana. The network was based on statistical parsimony for the 60 haplotypes retrieved from the COI sequences. Numbers indicate haplotype identity. Colors on network correspond to each locality as shown in the legend.

Genetic variation descriptors (Table 2) show that within populations the number of segregated sites range from 1 to 16. The highest nucleotide variation value was found in the population of Neblinas (π = 0.0149, S = 7). This locality had nearly 30 times the nucleotide variation found in Flor de Chiapas (π = 0.0005, S = 1). The lowest value for haplotype diversity was at Lachiguiri (Hd = 0.2), while the highest value (i.e., 1) was at Acahuizotla (Guerrero), Neblinas (Queretaro) and Tepexco (Puebla). The global Fu’s F value (i.e., -7.935; p < 0.001) is consistent with a significant demographic expansion history. However, none of the individual populations showed statistically significant values for this descriptor.

The SAMOVA (k change from 2 to 15) showed a gradual increase in the F_CT values (Table 3). The highest F_CT was reached with 15 groups, where the genetic variability among groups reached approximately 80%. This result suggests that practically every locality is an isolated group. The only two groups recovered were: Xalapa-La Orduña-San Antonio Paso del Toro (central Veracruz state) and Acahuizotla-Tepexco. The Mantel test revealed a significant positive correlation (r = 0.3756, p < 0.0001), which suggests isolation by distance. The AMOVA showed also significant genetic structure (F_ST = 0.77, P < 0.0001). Approximately 77% of molecular variance occurs among populations, while the remaining 23% is explained by the genetic differences within populations (Table 4). Population pairwise F_ST values revealed significant differentiation between most populations (Table 5). Gene flow estimations among populations were extremely varied, fluctuating from M = 0.000 (i.e., between Flor de Chiapas-San Antonio Paso del Toro) to a continuous exchange of individuals between Xalapa-La Orduña, and Xalapa-Cerro Colorado.

The best-fitting nucleotide substitution model was TrN+G+I. The model’s specific parameters were: base frequency = 0.3175, 0.1536 0.1489, 0.38; nst = 2; ts:tv ratio = 4.7931; Gamma shape = 0.5600; ncat = 4; and pinvar = 0.7490. The trees recovered by Bayesian (Fig. 3) and ML (See Suppl. material 4) inferences have congruent topologies, with slight changes in the internal relationships. The base of the F. mexicana clade is a polytomy formed by individuals of populations mainly distributed in the Sierra Madre del Sur, Pacific and Trans-Mexican Volcanic Belt, a Flor de Chiapas clade (FDC [PP = 1]), and a “big clade” (PP = 1) that encompasses the populations from Tierras Altas de Chiapas (TAC) to Sierra Madre Oriental (SMO). This “big clade” presents a substructure recovering a monophyletic group formed by Valles Centrales de Oaxaca (SMS: VCO) which is the sister group to the clade containing the rest of the sampled populations.

Figure 3. 

Genealogy recovered based on Bayesian inference. Numbers indicate the Bayesian posterior probabilities. The four main lineages are Flor de Chiapas (FDC), Sierra Madre del Sur that encompasses the Valles Centrales de Oaxaca (SMS: VCO), Tierras Altas de Chiapas (TAC), and Sierra Madre Oriental (SMO).

The chronogram (Fig. 4) showed a radiation event during the Middle Pleistocene. The results showed that F. mexicana split from the outgroup approximately 450,000 years ago. The diversification events that originated the main lineages within the species occurred nearly 150,000 years ago. The skyline-plots showed that the populations were kept in stasis for the most part of the Pleistocene, followed by a demographic expansion period towards the end of this epoch. For the lineage of TAC, a small increase in Ne occurred 3,500 years ago. In the SMO lineage, the increase of Ne occurred approximately 1,000 years ago. For the species (global analysis) the skyline-plot showed a demographic increment process that started approximately 4,000 years ago (Fig. 4).

Figure 4. 

Chronogram and Skyline-plots based on COI sequences of Falagonia mexicana. For the chronogram the time is in millions of years. The Skyline-plots are shown for the clades going through demographic changes (SMO: Sierra Madre Oriental, TAC: Tierras Altas de Chiapas, and Global: Whole genealogy). Skyline-plots show mean Ne, plus 95% HDP confidence limit.

We thank Dr. Domínguez Martínez and Dr. Sánchez Palmeros for their pertinent comments on this paper, as well as Jonathan Sulvarán-Salazar for his assistance in field. QJSJ is indebted to Enio Cano and Eunice E. Echeverria for the facilities provided to collect in Guatemala and El Salvador, respectively. We are grateful to A. Valencia for helping us with the English version of this paper. Dr. Janet Nolasco-Soto provided technical and logistical aid during the lab work. We thank Dr. J. L. Navarrete-Heredia and Dr. L. Delgado for the donation of biological material from La Concordia and Tepexco, respectively. In addition, we thank Dr. Gusarov for allowing us to use the outgroup’s sequence. Similarly, we thank the anonymous reviewers for their comments, as well as all the people who in any way helped us during the development of this research.

This work was partially supported by a CONACyT (Consejo Nacional de Ciencia y Tecnología) fellowship (# 164479) to QJSJ. Complementary financial support was obtained by the program “Apoyos Económicos para Estudiantes de Posgrado” (Instituto de Ecología, A. C.). QJSJ is grateful with IdeaWild organization for providing some equipment used during fieldwork and thanks the Universidad Veracruzana for the support to carry out a sabbatical stay at the Instituto de Ecología, A. C. to finish this manuscript. The funders were not involved in anyway regarding the study design, data collection and analysis, decision to publish, or preparation of the manuscript.