Co1 DNA supports conspecificity of Geomyphilus pierai and G. barrerai (Coleoptera, Scarabaeidae, Aphodiinae) and is a good marker for their phylogeographic investigation in Mexican mountains

Abstract Members of Geomyphilus are associated with rodent burrows, such as pocket gophers and prairie dogs. In Mexico, they are found in the mountains of the Mexican Volcanic Belt and in Sierra Madre Oriental. Our study aims to initiate the exploration of the dispersal modes of Geomyphilus pierai and Geomyphilus barrerai from burrows of pocket gophers. In order to estimate the dispersal scale of the beetles, the utility of mitochondrial and nuclear molecular markers for studying the phylogeographic structure of this complex of species (Geomyphilus pierai and Geomyphilus barrerai) was tested from 49 beetle individuals. High intraspecific and intra-mountain nucleotidic diversity was captured from this sample using Co1 mitochondrial sequences, whilst the ITS2 nuclear ribosomal sequence did not allow observing informative variation. Mitochondrial phylogenetic analysis revealed that the specific delineation between the two species under study was doubtful. In this preliminary study, Co1 was shown to be a good marker for elucidating dispersal routes of the burrowing rodent-associated beetles.


Introduction
The subfamily Aphodiinae is one of the most diversified groups within the Scarabaeoidea (classification sensu Dellacasa et al. 2001) showing a worldwide distribution (Dellacasa 1987(Dellacasa , 1988a(Dellacasa , b, 1991(Dellacasa , 1995 and having a remarkably high generic and specific diversity. The Aphodiinae are the predominant group in the communities of the dung beetles in the cold temperate areas of the Palearctic region (Hanski 1991) and also in the Nearctic region (Lobo 2000).
Members of Geomyphilus are endemic of the Nearctic region and are associated with rodent burrows. Common hosts include pocket gophers (Geomyidae) and prairie dogs (Sciuridae: Cynomis), but some data indicate that rarer Aphodiinae species may be found in association with kangaroo rats (Heteromyidae: Dipodomys), voles (Cricetidae: Microtus), or other rodents . In Mexico, members of Geomyphilus are found in the mountains of the Mexican Volcanic Belt (MVB) and in Sierra Madre Oriental, including the two presumed sibling species: Geomyphilus pierai (Lobo & Deloya, 1995) and G. barrerai (Islas, 1955).
While the habits and distribution patterns of their associated rodents are rather well-known (Demastes et al. 2002;Hafner et al. 2005;Fernández et al. 2014), these beetles' dispersal modes were unknown to date. Most of the southern pocket gophers, such as Cratogeomys merriami and Thomomys umbrinus, live at high elevations in the Mexican mountains. Their burrowing activity results in a net of galleries under the ground in many habitats of these mountains. The two species we are focusing on (G. pierai and G. barrerai) are known to live in close association with pocket gophers, with the beetles recurrently found nesting in and feeding on the rodent feces Halffter 1994, Deloya andLobo 1995, Arriaga-Jiménez pers. obs.). These nidicolous species were up to now only found underground, inside the rodent nests, latrines or galleries, with very few exceptions. Geomyphilus beetles possess a priori functional wings, which could allow a long distance spread by flying. However, the fact that the beetles live in a close association with rodents and that they were mainly found underground suggests that these scarabs may present mainly a short distance spread by going through the gophers' galleries. Alternatively, the rarity of the collecting outside the burrows might be a result of undersampling due to the short flight period of the beetles.
We undertook a preliminary molecular study to investigate Geomyphilus dispersal modes and clarify whether the preferred route is underground or aerial. Measuring the amount of gene flow between the populations in different mountains could help understanding the spatial scale of the Geomyphilus dispersal. Indeed, one may presume that if the species are mainly moving within the rodent galleries, the gene flow between the populations in different mountains should be less than if they fly outside.
To initiate the exploration of the dispersal routes of coprophagous scarabaeids narrowly associated to rodents, the development of molecular genotyping is needed, allowing to characterize the population genetic structure at different spatial scales within at least one species living in rodent burrows. Candidate molecular markers should provide genetic variation detectable at spatial scales appropriate to assess aerial dispersal. We assumed that the aerial dispersal might allow individuals to quickly fly from one mountain to another, whilst the underground dispersal most likely would limit the beetle circulation to a single mountain or at most to the interconnected mountains. Therefore, we decided to work on the individuals sampled from the nests located in the different mountains of MVB.
Moreover, several species of the genus Geomyphilus can be found together in the same gopher nest (Deloya and Lobo 1995, Arriaga-Jiménez, pers. obs.) and the two above cited species were recurrently found together in rodent nests by AAJ (during her PhD fieldwork). Because of such a co-occurrence and due to the rather difficult morphological discrimination between them, the interspecific molecular delineation was needed in first.
Here we present the results from a preliminary molecular investigation of G. pierai and G. barrerai aiming at (1) assessing interspecific molecular divergence and (2) testing molecular markers to allow phylogeographic studies to be performed in these species. We sequenced one mitochondrial and one nuclear DNA regions from individual beetles sampled from the burrows of pocket gophers (C. merriami and T. umbrinus) in four mountains of the MVB. The mountains under study were: La Malinche, Cofre de Perote, Pico de Orizaba and Sierra Negra (Fig. 1). The first one is separated from the others by the High Plateau, while the other three have some connectivity, especially Pico de Orizaba and Sierra Negra which are located next to each other.

Methods
To collect Geomyphilus beetles, the gopher nests were excavated (only those that seemed recently used). The galleries were followed up to find the nests and latrines where the dung beetles might be found. This sampling technique is very effort-consuming because the nests can be located very deep and they are not always found. Once the nests and latrines were reached, the organic matter was removed as well as the soil within 20 cm from it, and then all the soil was carefully inspected for scarabs. The depth of the nest and latrines varies depending on the type of soil, slope and surrounding vegetation. They were found from 30 cm to 1.5 m deep in different biotopes (scrublands and grasslands) in the studied mountains (Fig. 2). All the individuals were immediately transferred to 90% ethanol and conserved for further morphological identification and molecular analyses. All specimens examined were identified by Giovanni and Marco Dellacasa.
The DNA was extracted from ethanol-preserved specimens. Either of two different techniques were used, both designed to obtain molecular information while keeping voucher cuticle of each individual. For the first one, after complete evaporation of ethanol, the cuticle of each specimen was cut at two points on the prothoracic-mesothoracic joint soft integument in order to facilitate the lysis of internal tissues. This was done with a sterile syringe needle on a glass slide before insertion of individual in a vial for proteinase K digestion. For the second technique, one leg was removed from the beetle, then put in an Eppendorf vial and crushed with a pestle. For both techniques, the entire cut individual or the crushed leg were immersed in 100 ml Buffer ATL + 15 ml Proteinase K and incubated for 19-30 h at 70 °C, before proceeding to the DNA purification using Qiagen's DNeasy Blood & Tissue Kit following the manufacturer's protocol. The specimens' cuticles (tech. 1) were conserved in glycerol and the body remaining (tech. 2) in ethanol as vouchers.
Checking congruence between phylogenetic topologies from genes with different transmission patterns such as mt DNA (maternal transmission) and nuclear DNA (biparental transmission) is useful to confirm cryptic species (see Donelly et al. 2013). For this purpose, Co1 and IT2 appeared to be good candidates, since they were shown to be variable both at the inter-and intraspecific levels in some scarab beetle genera (Wirta 2009). Two fragments from one mitochondrial protein-coding gene (Co1) and one nuclear ribosomal internal transcribed spacer (ITS2) gene were amplified. Primer pairs used for amplification of Co1 were LCO_1490 and HCO_2198 (Folmer et al. 1994) or LEP-F1 and LEP-R1 (Hebert et al. 2004); and for the nuclear region ITS2 we used ITS2f (Navajas et al. 1998) and RhITSR (De Rojas et al. 2002).
Polymerase chain reaction (PCR) was performed using the Qiagen's Taq PCR Core kit. The following reagent concentrations were used for a final volume of 25 µL per tube: 1×PCR buffer, 0.036U/µl Qiagen Taq Polymerase, 300 µM dNTPs and 0.6 µM of each Co1 primer or 0.43 µM of each ITS2 primer, and 2.5 mM MgCl 2 for CO1 or 3.0 mM MgCl 2 for ITS2. For the Co1, thermal cycling parameters comprised an initial denaturation step of 10 minutes at 94 °C, followed by 5 cycles of 94 °C for 40 s, 49 °C for 60 s and 62 °C for 60 s, and 35 cycles of 94 °C for 40 s, 52 °C for 60 s and 62 °C for 60 s, with a final elongation at 62 °C for 10 min. For the ITS region, thermal cycling parameters comprised an initial denaturation of 10 minutes at 94 °C, followed by 5 cycles of 94 °C for 40 s, 56 °C for 60 s and 72 °C for 60 s, and 35 cycles of 94 °C for 52 s, 52 °C for 60 s and 72 °C for 60 s, with a final elongation at 72 °C for 10 min. Purification and sequencing was realized by Genoscreen (France, Lille) using a 96-capillary sequencer ABI3730XL.
Obtained DNA sequences were aligned using Muscle (Edgar 2004) in Seaview4.1 (Galtier et al. 1996). The absence of any stop codon was first checked in the Co1 sequence alignment. Maximum likelihood phylogenetic analyses were performed using PhyML (Guindon and Gascuel 2003) and nucleotidic diversities were compared between samples and between taxa using DNAsp 5.10.01 (Rozas and Rozas 1995;Librado and Rozas 2009). A median-joining network (Bandelt et al. 1999) of DNA fragments was constructed using NETWORK v4.5.1.6 (fluxus-engineering.com) to visualize the frequency and geographical distribution of Co1 haplotypes.

Results and discussion
In order to estimate the utility of the two DNA markers to assess the phylogeographic structure of this complex of species (G. pierai and G. barrerai; Fig. 3), a 700-bp Co1 fragment was obtained from 28 individuals and a ca. 830-bp ITS2 fragment from 47 individuals collected in 21 rodent nests distributed across the four different mountains of the MVB (total 49 sequenced individuals, see Online supplementary materials). Intraspecific and intra-mountain diversity was captured from this sample using mitochondrial sequences, not nuclear sequences (Table 1). Indeed, 17 Co1 haplotypes were recorded among individuals identified as belonging to either G. pierai or G. barrerai, with 46 polymorphic sites, a 4.1 mean pairwise percent of nucleotidic differences and a 0,963 Hd value (haplotype diversity). In contrast, 4 ITS2 haplotypes were recorded, with only 3 polymorphic sites ( Figure  4). Not only differences were very scarce among the nuclear DNA fragments, but also a single haplotype was largely dominating, contrarily to what was observed in some other scarab species (e.g. Wirta 2009). No clear association between any mitochondrial lineages was noticed and the largely dominant haplotype was also found in the three outgroup individuals (N. perotensis). As a result, we considered this marker inappropriate to investigate phylogeographics in Geomyphilus spp. and only performed subsequent analyses on the Co1 dataset.
Species boundaries between G. barrerai and G. pierai were not confirmed either by mitochondrial or nuclear data. When considering concordance between morphological taxonomic assignation and molecular information, not only none of the 46 Co1 polymorphic sites and none of the 4 ITS2 polymorphic sites did segregate accordingly, but also, the distribution of the two taxa across the Co1 topology was dispersed and intermingled, without any clade grouping together all individuals of any of them (Fig. 4). Although some other cryptic species (different than either G. barrerai or G. pierai) could eventually be revealed by using a more rapidly-evolving nuclear gene, the current specific delineation is most likely unable to match any actual sexual isolation between these current taxa. If two or more other cryptic lineages were confirmed by some other nuclear analyses, a wide morphological screening involving genotyped individuals from diverse origins would be required to improve specific description and allow determining whether any characters are both discriminant between these entities and stable within each of them. This is strengthened by the fact that higher nucleotidic diversity was found in Co1 sequences between mountains than between taxa (see Table 2). All that  strongly suggests that these two taxa represent either a single biological species or a differently structured species complex. As a result, G. pierai should be considered for synonymization with G. barrerai following the ICZN principle of priority. Subsequently considering the sample as monospecific, so much variation in the Co1 sequence dataset, obtained from the present locally-restricted sample is promising, providing basic tools for further analyses. Besides, several different mitochondrial lineages were found within the beetle sample, with no clear correspondence between mountains and genetic structure (Fig. 5). Although a group of four close haplotypes could appear to be if not mountain specific, more or less place specific (see the right orange and brown clade sampled from the two closest mountains Cofre de Orizaba and Sierra Negra in Fig. 5), all the main clades group together haplotypes from two or three more or less distant mountains. More importantly, some Co1 haplotypes were found to be shared between individuals sampled from different mountains, including a priori non burrow-connected mountains (Malinche and Pico de Orizaba). Such an apparent lack of geographical structure at the Mexican Volcano Belt scale in a so diverse sequence dataset may suggest either that aerial spread commonly occurs or that some cryptic sub-structure exists, though no conclusion may be drawn from a so small sample. These preliminary results provided a first insight of the phylogeographic structure of the G. pierai/G. barrerai complex. A further study using at least Co1 as a marker on a much larger sample should help to elucidate dispersal routes and gene flow of these nidicolous rodent-associated beetles.