The role of chromosomal rearrangements and geographical barriers in the divergence of lineages in a South American subterranean rodent (Rodentia: Ctenomyidae: Ctenomys minutus)

Abstract

Identifying factors and the extent of their roles in the differentiation of populations is of great importance for understanding the evolutionary process in which a species is involved. Ctenomys minutus is a highly karyotype–polymorphic subterranean rodent, with diploid numbers ranging from 42 to 50 and autosomal arm numbers (ANs) ranging from 68 to 80, comprising a total of 45 karyotypes described so far. This species inhabits the southern Brazilian coastal plain, which has a complex geological history, with several potential geographical barriers acting on different time scales. We assessed the geographical genetic structure of C. minutus, examining 340 individuals over the entire distributional range and using information from chromosomal rearrangements, mitochondrial DNA (mtDNA) sequences and 14 microsatellite loci. The mtDNA results revealed seven main haplogroups, with the most recent common ancestors dating from the Pleistocene, whereas clustering methods defined 12 populations. Some boundaries of mtDNA haplogroups and population clusters can be associated with potential geographical barriers to gene flow. The isolation-by-distance pattern also has an important role in fine-scale genetic differentiation, which is strengthened by the narrowness of the coastal plain and by common features of subterranean rodents (that is, small fragmented populations and low dispersal rates), which limit gene flow among populations. A step-by-step mechanism of chromosomal evolution can be suggested for this species, mainly associated with the metapopulation structure, genetic drift and the geographical features of the southern Brazilian coastal plain. However, chromosomal variations have no or very little role in the diversification of C. minutus populations.

Introduction

The role of chromosomal rearrangements in the speciation process has been widely discussed over several decades. Traditional chromosomal speciation models propose that the fitness of heterozygous carriers is diminished owing to meiotic abnormalities, leading to reductions in gene flow between chromosomally divergent populations and thus facilitating reproductive isolation (White, 1978). However, recent models suggest that chromosomal changes reduce gene flow and accelerate the divergence between populations by reducing recombination in heterokaryotypes and extending the effects of linking isolated genes or locally adapted alleles, and not because of underdominance (Rieseberg, 2001). However, it is questionable whether all the kinds of chromosomal rearrangements could be involved in the speciation process (Baker and Bickham, 1986).

On the other hand, the role of geographical factors in the process of biological diversification is increasingly understood, and its contribution as an evolutionary force shaping the genetic structure of the species is broadly accepted (Barraclough and Vogler, 2000). The geographical distribution of a species depends on, among other intrinsic and extrinsic factors, when and where it originated and what kinds of geographical barriers it encountered. Mountain chains, deserts, habitat fragmentation or watercourses can sometimes be insuperable geographical barriers for individuals that reach them, reducing or preventing gene flow between populations located on opposite sides, and thus accelerating the process of divergence (Slatkin, 1987). Steinberg and Patton (2000) suggested that the geographical context on the process of speciation in subterranean rodents is trivial, in which the vicariance or peripheral isolation may account for the allopatric divergence.

The burrowing rodents of the genus Ctenomys, popularly known as tuco-tucos, comprise approximately 60 living species, showing impressive intra- and interspecific karyotypic variation (Reig et al., 1990; Woods and Kilpatrick, 2005). Some of their characteristics, such as low rates of adult dispersal (Busch et al., 2000) and their distribution in relatively small and fragmented populations, foster the establishment of small genetic units where genetic variation is low and interpopulation divergence is high (Wlasiuk et al., 2003). These characteristics make the species more vulnerable to stochastic demographic factors and changes in their habitats (Patton et al., 1996), and also favor the fixation of new chromosomal rearrangements, which have been considered to have an important role in the diversification of ctenomyids (Reig et al., 1990; Lessa and Cook, 1998). This genus provides an excellent study system for exploring the role of discontinuities in the environment and chromosomal rearrangements in the speciation process (Mora et al., 2006, 2007; Fernández-Stolz, 2007; Tomasco and Lessa, 2007) and ultimately for identifying the important genes involved in reproductive isolation of these taxa. Thus, the investigation of population genetic structure presented in this manuscript is a first step in the process of exploring these large questions.

Ctenomys minutus Nehring, 1887 has a narrow distribution on the southern Brazilian coastal plain. The landscape of the region is a mosaic of innumerable lakes and lagoons, rivers, swamps, sand fields and dunes. During the Quaternary Period, fluctuations in the Atlantic Ocean sea level produced large lateral displacements of the shoreline, originating four barrier–lagoon systems that shaped the present coastline. The first three barrier–lagoon systems originated during the Pleistocene Period (hereafter called Barriers I, II and III), and the fourth system was formed in the Holocene Period (Barrier IV; Supplementary Figure S1; Tomazelli et al., 2000).

C. minutus inhabits only the sand fields (corresponding to Barriers II and III) in the southern portion of its distribution; and in the north, these tuco-tucos preferentially occupy the first dune line (Barrier IV; Figure 1 and Supplementary Figure S1; Freygang et al., 2004). The range includes many watercourses that have existed in different periods, affecting the local rodent populations over different time scales. According to studies conducted by Weschenfelder et al. (2008) and Baitelli (2012), three major paleochannels existed on the coastal plain (Paleojacuí north, Paleojacuí south and Paleocamaquã), but they were completely closed by sediments during the Pleistocene and Holocene Periods. Nowadays, the Patos Lagoon and the Araranguá and Mampituba Rivers form major discontinuities in the landscape of the region, and all of them flow into the Atlantic Ocean (Figure 1).

Figure 1

Sampling localities of Ctenomys minutus, and the distribution of the parental karyotypes and intraspecific hybrid zones. The seven mtDNA haplogroups are highlighted in the figure. BAR, Barco Beach; BJ1, Bujuru 1; BJ2, Bujuru 2; EBL, East Barros Lake; FOR, Fortaleza Lake; FSM, Farol de Santa Marta; GAI, Gaivota Beach; GUA, Guarita Beach; ILH, Ilhas; JG, Jaguaruna; MC, Morro dos Conventos; MOS, Mostardas; PAL, Palmares do Sul; PAS, Passinhos; PIT, Pitangueira; PST, Passo de Torres; OSO, Osório; 35, road km 35; 53, road km 53; 64, road km 64; 96, road km 96; 108, road km 108; 115, road km 115; SBL, South Barros Lake; 17S, 17 km south of Mostardas; 26S, 26 km south of Mostardas; SJN, São José do Norte; TRA, Tramandaí; TV1, Tavares 1; TV2, Tavares 2.

C. minutus is a highly karyotype–polymorphic species, with diploid numbers ranging from 42 to 50 and autosomal ANs from 68 to 80. The chromosomal variation consist of eight parapatrically distributed parental karyotypes (2n=50a, 48c, 46a, 48a, 42, 46b, 48b and 50b) and six intraspecific hybrid zones, among them: (i) 50a × 48c=49a; (ii)46a × 48a=47a; (iii)48a × 42=43, 44, 45, 46; (iv) 42 × 46b=43, 44, 45; (v) 46b × 48b=47b; and (vi) 48b × 50b=49b. The intraspecific hybrid zones were determined based on diploid numbers and chromosomal rearrangements observed through G-band patterns compared with the adjacent parental karyotypes (Freitas, 1997; Gava and Freitas, 2002, 2003; Freygang et al., 2004; Castilho et al., 2012; this study). The distribution of these karyotypes shows 2n=50 at both ends of the geographical range, whereas toward the middle of the distribution, it is progressively reduced to 2n=42 through Robertsonian rearrangements, tandem fusions/fissions and a pericentric inversion, which originated the distinct karyotypic systems ‘a/c’ and ‘b’; see Figures 1 and 2 (Freitas, 1997; Freygang et al., 2004).

Figure 2

Comparison of chromosomal rearrangements found among parental karyotypes based on the G-band patterns. The karyotype 2n=50a, AN=76 was used as standard to compare the chromosomal rearrangements in the other karyotypes. Squares highlight the main chromosomal rearrangements involved in the intraspecific hybrid zones. Modified from Freitas (1997), Gava and Freitas (2003) and Freygang et al. (2004).

Considering the several potential water barriers acting on different time scales along the southern Brazilian coastal plain, and since tuco-tucos seem to be very poor swimmers (Reig et al., 1990), we investigated here if perennial streams could lead to a process of genetic differentiation between C. minutus populations in allopatry and also if geographical features could have influenced the fixation of chromosomal differences. We also explored the possibility that the chromosomal rearrangements found in this species could represent incipient genetic barriers to reduce gene flow between chromosomally divergent populations. In this contribution, we attempted to examine both historical and contemporary processes that shaped and continue to shape the geographical genetic structure of C. minutus, regarding the coastal plain dynamics, by means of the mitochondrial DNA (mtDNA) control region (CR) and cytochrome c oxidase subunit I sequences (COI), 14 microsatellite loci, using specimens sampled and karyotyped on a fine scale.

Materials and methods

Sample collection and cytogenetics

The samples analyzed included 3 to 25 specimens of C. minutus from each of 30 localities, totaling 276 individuals analyzed for mtDNA, and 340 individuals analyzed for microsatellite loci, across the species’ entire distributional range (see Table 1 and Figure 1). Most of the specimens were collected and karyotyped previously in studies by Freitas (1997), Gava and Freitas (2002, 2003, Freygang et al. (2004) and Castilho et al. (2012) (Table 1). The tissue samples were preserved in 95% ethanol and stored at −20 °C in the collection of the Laboratório de Citogenética e Evolução of the Departamento de Genética of the Universidade Federal do Rio Grande do Sul, Brazil. Additional specimens were collected and karyotyped for this study and are highlighted in Table 1. For these specimens, the mitotic preparations were obtained from bone marrow, according to Ford and Hamerton (1956). The diploid and autosomal numbers were determined by analyses of at least 20 metaphase spread cells stained with Giemsa.

Table 1 Genetic structure and variability observed among the localities (abbreviations within parentheses) sampled for Ctenomys minutus for both mitochondrial DNA sequences and nuclear microsatellite markers

DNA amplification, sequencing and genotyping

Total DNA was extracted using a standard phenol:chloroform protocol (Sambrook and Russel, 2001). Two fragments of mtDNA were amplified by polymerase chain reaction: Part of HVS1 from the CR was amplified using the primers TucoPro (5′-TTCTAATTAAACTATTTCTTG-3′; Tomasco and Lessa, 2007) and TDKD (5′-CCTGAAGTAGGAACCAGATG-3′; Kocher et al., 1989), following a protocol modified from Tomasco and Lessa (2007). The cytochrome c oxidase subunit I (COI) followed the protocols suggested in http://www.barcoding.si.edu/DNABarCoding.htm, using the primers LCO-1490 (5′-GGTCAACAAATCATAAAGATATTGG-3′) and HCO-2198 (5′-TAAACTTCAGGGTGACCAAAAAATCA-3′; Folmer et al., 1994).

The polymerase chain reaction products were purified using Exonuclease I and Shrimp Alkaline Phosphatase (GIBCO-BRL Life Sciences/Invitrogen, Carlsbad, CA, USA), following the guidelines of the suppliers, and sequenced in an ABI 3730 (Applied Biosystems, Foster City, CA, USA) automated sequencer, using the forward primers TucoPro and LCO-1490.

In addition, 14 microsatellite loci were amplified, using fluorescently labeled primers, isolated from the Argentinean species C. haigi (Hai2, Hai3, Hai4, Hai5, Hai6, Hai9, Hai10, Hai12; Lacey et al., 1999) and C. sociabilis (Soc1, Soc2, Soc3, Soc4, Soc5, Soc6; Lacey, 2001). Polymerase chain reaction amplifications were carried out following the protocols described by Lacey et al. (1999) and Lacey (2001). The genotypes were obtained in the sequencer ABI 3730, mixing two primer pairs with different labels, per reaction. To define the allele sizes, we used the program PeakScanner 1.0 (http://www.appliedbiosystems.com).

Data analyses

Genetic diversity

Sequence electropherograms were visually inspected using Chromas 2.33 (http://www.technelysium.com.au/chromas.html), and aligned using the Clustal W algorithm with default options, implemented in Mega 4.0.2 (Tamura et al., 2007). All analyses were performed with the three data sets: CR and COI separately, and a concatenated data set of CR+COI (CC); however, mainly the CC results will be presented here. Measures of mtDNA diversity, such as average number of nucleotide differences (k), average number of nucleotide differences per site between all pairs of sequences (π), definitions of haplotypes (H) and haplotype diversity (H_d), were calculated with DnaSP 5.00.03 (Librado and Rozas, 2009).

We used Arlequin 3.5.1.2 (Excoffier and Schneider, 2005) and Genepop 4.0 (Rousset, 2008) to test deviations from the Hardy–Weinberg equilibrium, to calculate the observed and expected heterozygosis (H_o and H_e) and to perform the analyses of linkage equilibrium (LE) for microsatellite loci. The significance levels (α=0.05) of Hardy–Weinberg equilibrium and LE were adjusted using sequential Bonferroni corrections. Genotyping errors were tested with Micro-Checker (van Oosterhout et al., 2004), under a 95% confidence interval. The measurements of microsatellite diversity were assessed using Arlequin 3.5.1.2.

Genetic differentiation

Wright’s F-statistics were applied to analyze the within- (F_IS) and between- (F_ST) population structures for microsatellite data, using Arlequin. Pairwise R_ST population estimates (Slatkin, 1995) was also calculated, using Genepop 4.0. Localities with fewer than five specimens were excluded from pairwise comparisons, and sequential Bonferroni corrections were applied to adjust the significance levels (α=0.05).

The program Arlequin 3.5.1.2 was used to perform analyses of molecular variance (AMOVA) to investigate two different scenarios of hypothesized population subdivisions, regarding the karyotypic variation and the potential geographical barriers. All the specimens collected in intraspecific karyotypic hybrid zones were excluded from both analyses. First, the samples were subdivided into eight groups according to parental karyotypes: (i) only specimens with 2n=50a, from Farol de Santa Marta (FSM) (the complete names of the localities and the corresponding abbreviations are listed in Table 1); (ii) specimens with 2n=48c, from Ilhas (ILH); (iii) 2n=46a, all specimens from Morro dos Conventos (MC), Gaivota Beach (GAI), Passo de Torres (PST), Guarita Beach (GUA), Barco Beach (BAR), Tramandaí (TRA) and Osório (OSO); (iv) 2n=48a, from Passinhos (PAS), Pitangueira (PIT), Palmares do Sul (PAL), road km 35 (35), road km 53 (53), road km 64 (64) and road km 96 (96); (v) 2n=42, from road km 115 (115) and Mostardas (MOS); (vi) 2n=46b, from Tavares 1 (TV1); (vii) 2n=48b, from Bujuru 1 (BJ1); and (viii) specimens with 2n=50b, from São José do Norte (SJN). The same data set was used in a second test, but the samples were separated into seven groups, considering present and possible former geographical barriers: (i) specimens from FSM and ILH, isolated to the south by the Araranguá River; (ii) specimens from MC, GAI and PST, isolated to the south by the Mampituba River; (iii) specimens collected in GUA and BAR, delimited to the south by the transition between sand fields and dunes; (iv) all specimens from TRA, OSO, PAS, PIT, PAL and 35, isolated in the past to the south by the Paleojacuí north; (v) samples from 53, 64, 96, 115 and MOS, isolated in the past to the south by the Paleojacuí south; (vi) TV1, isolated in the past to the south by the Paleocamaquã; and (vii) specimens sampled in BJ1 and SJN, isolated in the southernmost part of the C. minutus distribution.

The F_ST and R_ST estimates, and the AMOVA tests were also performed using reduced data sets, excluding loci showing evidence of potential null alleles and/or linkage disequilibrium, to ensure that these loci do not bias the estimations of gene flow and the genetic structure results observed for C. minutus.

Isolation-by-distance analysis

To search for positive correlations between genetic and geographical distances, a simple Mantel test was implemented, based on pairwise estimates of F_ST and linear geographical distances between localities. Partial Mantel tests were also computed to observe whether genetic and chromosomal variations are positively correlated, while controlling for the spatial influence. These analyses were computed in Arlequin 3.5.1.2.

The Robertsonian rearrangements and the pericentric inversion shown in Figure 2 were coded as shown in Supplementary Table S1. As we do not have the G-band patterns of some specimens from locality road km 108 (108), and all specimens from 17 km south of Mostardas (17S) and 26 km south of Mostardas (26S), we excluded these individuals from the analyses. The chromosome rearrangements indicated for karyotype 48c were inferred based on its diploid and autosomal arm numbers (ANs), and also on the G-band patterns of karyotypes 50a and heterozygotes 49a. The pairwise karyotypic distances between sampling sites were calculated based on Nei’s genetic distance (Nei, 1978), implemented in Popdist 1.2.4 (Guldbrandtsen et al., 2000).

Moreover, a spatial autocorrelation analysis was performed in the program Alleles In Space 1.0 (Miller, 2005), assuming 30 equal distance classes with unequal sample sizes. We also checked whether the chromosomal rearrangements follow a pattern of isolation by distance. The matrix of chromosomal rearrangements was applied in the program Alleles In Space as codominant diploid data, together with the geographical coordinates of the specimens, to implement the spatial autocorrelation analysis and Mantel test analyses. The statistical significances of the tests were assessed assuming 1000 random permutations

Phylogenetic analysis

An incongruence length difference test was implemented in the program PAUP* 4.0 (Swofford, 2002), to ensure that both CR and COI data do not have significantly different phylogenetic signals, and can be combined. We applied a heuristic search, under 1000 random permutations. We assumed a significance level of 0.01, as the incongruence length difference test is prone to the type I error (Cunningham, 1997). A phylogenetic analysis using Bayesian inference was performed with MrBayes 3.1.2 (Huelsenbeck and Ronquist, 2001). The CR and COI sequences were combined into a single data set (incongruence length difference, P=0.015), with the data partitioned by fragment, for which the appropriate model of nucleotide sequence evolution, determined in MrModeltest 2.2 (http://www.abc.se/~nylander), was applied (CR: HKY+I+G; and COI: HKY+G). Two independent replicates of the Markov chain Monte Carlo search were performed with 1 000 000 simulations, each with four chains, sampling trees every 100 generations. After the convergence test, the first 500 trees were discarded. Homologous sequences of C. torquatus Lichtenstein, 1830 and C. flamarioni Travi, 1981 were used as the outgroup. The topological relationship between the haplotypes was estimated using the program Network 4.5.1.0 (http://www.fluxux-engineering.com) with the median-joining approach, for the concatenated mtDNA data set.

To better visualize the relationships among parental karyotypes represented by a dendrogram, we implemented the neighbor-joining method in the program PAUP* 4.0, based on a matrix of presence/absence of chromosomal rearrangements, assuming 1000 random replicates.

Demographic history and divergence time estimates

The program Arlequin 3.5.1.2 was used to assess the demographic history of population expansion, employing Tajima’s D and Fu’s F_S neutrality tests, and mismatch distribution analysis. To test the goodness of fit of the observed mismatch distribution, we used the sum of squared deviations (SSD) and the Harpending’s raggedness index (RI) (Harpending, 1994). In addition, to examine past population dynamics and the time of the most recent common ancestor (tMRCA) for the seven main haplogroups identified in phylogenetic analysis, we applied the Bayesian skyline approach implemented in the program Beast 1.7.2 (Drummond and Rambaut, 2007), using the molecular models of evolution provided by MrModeltest 2.2. We applied a strict molecular clock with the substitution rates of 0.0213 (95% confidence interval: 0.0142–0.0287) and 0.0296 (95% confidence interval: 0.0165–0.0442) substitutions per site per million years for COI and CR, respectively, estimated for the genus Ctenomys by Roratto (2012) as normal priors, with 50 000 000 of Markov chain Monte Carlo simulations sampled every 1000 chains. Runs were summarized in TreeAnnotator 1.5.0 (Drummond and Rambaut, 2007), with a burn-in of the first 10% of the trees.

Bayesian clustering of populations

To determine the pattern of population genetic structuring and to assign individuals to populations, we used the Bayesian Markov chain Monte Carlo approach with and without spatial data, implemented in the programs Tess 2.3.1 (François et al., 2006) and Structure 2.3.3 (Pritchard et al., 2000), respectively. The number of clusters (K) was estimated in Tess 2.3.1 using the CAR admixture model, a burn-in of 50 000 sweeps followed by 50 000 iterations. We applied 10 independent runs for each K, ranging from 2 to 30. The values of deviance information criterion (DIC) of each run were plotted against the K_max, and the most likely number of clusters was considered when the DIC curve approaches a plateau. The Structure 2.3.3 analyses also followed 10 independent runs from K=2–30, 2 000 000 Markov chain Monte Carlo iterations, and a burn-in period of 1 500 000 steps, assuming an admixture model, correlated allele frequencies and a priori location information for each sample. The most likely number of clusters was based on the second-order rate of change in the log probability of the data (ΔK; Evanno et al., 2005).

We performed 50 additional runs for the most likely K found previously in Tess 2.3.1 and Structure 2.3.3, following the same parameters described above. We showed the results for the run that retrieved the highest mean of the estimated logarithm of probability of the data (ln Pr (X/K)) for Structure 2.3.3, and the lowest DIC value for Tess 2.3.1. We assumed an arbitrary threshold of membership probability q⩾0.80, to ensure that specimens with more than 80% of their genomic content represented in one cluster will be assigned to it.

Results

Cytogenetics

Of the 21 specimens collected at sampling sites 17S and 26S, 12 individuals were successfully karyotyped, revealing three specimens with a known karyotype (2n=42, AN=74), described by Freitas (1997), and nine specimens with four new karyotypes, intermediate between the forms 2n=42 and 46b, described in Supplementary Table S2 and Supplementary Figure S2A–D.

All specimens recovered in FSM had the karyotype 2n=50a, AN=76, previously described by Freitas (1997); the specimens from ILH had unknown karyotype, 2n=48c, AN=76 (Supplementary Table S2 and Supplementary Figure S2E); and at locality JG, we found specimens with 2n=48c and 50a, and heterozygotes with 2n=49a, AN=76.

Considering all cytogenetic data available for C. minutus, there is evidence that this species has nine diploid numbers (2n=42, 43, 44, 45, 46, 47, 48, 49, and 50), eight parental karyotypes (2n=42, 46a, 46b, 48a, 48b, 48c, 50a, and 50b) and six hybrid zones (between: 50a × 48c; 46a × 48a; 48a × 42; 42 × 46b; 46b × 48b; and 48b × 50b), comprising a total of 45 karyotypes known so far (Table 1).

Genetic diversity

From 398 bp of the CR, 43 polymorphic sites were obtained, resulting in a total of 40 haplotypes. The COI sequences were 620 bp in length, resulting in 44 variable sites and 35 haplotypes, and the concatenated data comprised 1018 bp, generating a total of 52 haplotypes. The levels of diversity in the three data sets were mostly moderate to high compared with indices of mtDNA diversity for other ctenomyids (overall estimates, π=0.0100–0.0236; k=6.20–15.59; H_d=0.947–0.965; Table 2).

Table 2 Mitochondrial genetic variability estimates for the seven main haplogroups identified for C. minutus and the total variation

For microsatellite data, the total number of alleles per locus among the 14 loci analyzed ranged from 7 to 16, and from 1 to 12 per locus per locality. Sampling sites FSM, JG, ILH, BAR, EBL, 35 and 96 showed one monomorphic locus; TV2, BJ1 and BJ2 showed two monomorphic loci; and SJN had six monomorphic loci. The monomorphic loci differed among the localities. The mean numbers of alleles per locality are presented in Table 1.

Departures from Hardy–Weinberg equilibrium were found in six sampling sites: OSO (Hai12 and Soc3); EBL (Hai4); SBL (Soc5); 96 (Soc3); MOS (Hai6); and BJ2 (Hai4 and Soc6). These deviations could be explained by (i) null alleles detected with Micro-Checker; (ii) as a consequence of reduced effective population sizes, which together with low rates of dispersal can lead to deviations from random mating, as evidenced in estimates of significant local F_IS in ILH, OSO, EBL, 96 and MOS (Table 1); or even (iii) the departures from Hardy–Weinberg equilibrium detected at EBL, SBL and BJ2 may be a consequence of hybridization events, as three of six loci that were not in equilibrium were from these hybrid localities. Pairwise comparisons of LE within localities revealed significant results among five pairs of loci: Hai3 × Soc5 (OSO), Hai10 × Soc2 (JG, OSO, 108), Hai12 × Soc1 and Soc4 × Soc5 (GAI) and Soc3 × Soc6 (MOS). However, these loci were not excluded from the subsequent analysis because there were no significant deviations from LE, for the same loci, in other pairwise comparisons for C. minutus. The influence of potential null alleles and loci under linkage disequilibrium seems to be negligible regarding the F_ST and R_ST estimates, and the AMOVA tests performed with reduced data sets, without these loci (data not show).

Phylogeographic pattern of haplotype distributions

Seven main clades were highlighted in the Bayesian phylogenetic tree (Figure 3a), all of them showing a strong relationship between specimen clustering and the geographical distribution. This same pattern could also be observed in the seven haplogroups retrieved by the median-joining haplotype network (Figure 3b). The northern localities were subdivided into two separate clades: North 1 was formed by all specimens from MC and 14 specimens from JG (six specimens with 2n=50a, three with 2n=49a and five with 2n=48c); and North 2 included the specimens from FSM, ILH and three specimens from JG (two with 2n=49a and one with 2n=50a). The middle of the C. minutus distribution was subdivided into three clades, among them: Coast, formed by specimens inhabiting only the shoreline (GAI, PST, GUA and BAR); Barros Lake, comprising sampling sites TRA, OSO, EBL, SBL, PAS, FOR, PIT and PAL, and one specimen from 35; and Mostardas, represented by the specimens from 35, 53, 64, 96, 108, 115 and MOS. The Tavares clade was represented by specimens from 17S, 26S, TV1 and TV2. Finally, the South clade, comprising specimens from BJ1, BJ2 and SJN.

Figure 3

(a) Bayesian phylogenetic tree; and (b) median-joining haplotype network topology obtained for concatenated mtDNA data. The seven main genetic haplogroups are highlighted by squares. Shading represents diploid numbers of parental karyotypes and intraspecific hybrids, as indicated in the legend. The abbreviations of the localities and the number of specimens per sampling site (within parentheses) are shown in the tree, and correspond to those in Table 1. The node posterior probabilities are given on the branches. Small black circles represent extinct or unsampled haplotypes in the network topology, and circle areas are proportional to the haplotype frequencies. *Diploid numbers of hybrids from crossing between 2n=42 × 48a; **diploid numbers of hybrids from crossing between 2n=42 × 46b. The color reproduction of this figure is available on the Hereditary journal online.

The mean tMRCAs calculated for each haplogroup were as follows: North 1—34 719 years (95% highest posterior density (HPD): 1613–82 431 years); North 2—179 390 years (95% HPD: 66 141–310 040 years); Coast—152 430 years (95% HPD: 58 048–259 190 years); Barros—122 130 years (95% HPD: 49 026–212 930 years); Mostardas—162 350 years (95% HPD: 65 339–276 080 years); Tavares—109 440 years (95% HPD: 37 531–194 680 years); and South—33 032 years (95% HPD: 1424–78 114 years).

In general, most haplotypes were separated by several mutational steps, and were limited to single localities, that is, private alleles. The shared haplotypes were mainly among neighboring sampling sites (Figures 1 and 3a and Table 1). Concerning chromosomal polymorphisms, only specimens from parapatric karyotypes and their hybrid forms shared haplotypes in the three mtDNA data sets (Supplementary Table S3).

The dendrogram of chromosomal rearrangements demonstrated that the karyotypes distributed at the north and south extremes of the geographical distribution of C. minutus (that is, 2n=50a, 48c and 50b) are more similar to each other than with the karyotypes observed in the middle of the distribution (Supplementary Figure S3).

Demographic history and isolation by distance

Although the plots of mismatch distributions were not clearly unimodal (Supplementary Figure S4), the tests of goodness of fit indicated that the observed and expected distributions of pairwise differences among haplotypes did not differ significantly under a model of recent population expansion, for COI (SSD=0.002, P=0.77; RI=0.006, P=0.88) and CC data sets (SSD=0.004, P=0.28; RI=0.004, P=0.24). In contrast, the observed and expected mismatch distributions differed significantly for the CR data set (SSD=0.012, P=0.001; RI=0.016, P=0.001).

There was no statistical support to accept a recent history of population expansion for C. minutus, according to the results of neutrality tests. This hypothesis is reinforced by high average number of nucleotide differences among haplotypes (ranging from k=6.20 to 15.59; Table 2), the limited distribution of haplotypes, and also by the absence of a star-like topology in the median-joining network, which is expected under a scenario of recent population expansion. The overall results of Tajima’s D showed positive values for CR and CC data and a negative D for COI, although none of the values was significant (Table 2). For the global Fu’s F_S test, all values were negative but not significant. With respect to the neutrality test results for the seven haplogroups of mtDNA, only Barros Lake showed negative and significant Fu’s F_S values, for the CC data set (Table 2). The effective population sizes were mostly constant among all the seven haplogroups, according to the Bayesian skyline analyses, and only some minor sudden retractions were observed from the last 13 000 years (Figure 4).

Figure 4

Bayesian skyline plots of concatenated mtDNA data for the seven main haplogroups of C. minutus: (a) North 1; (b) North 2; (c) Coast; (d) Barros Lake; (e) Mostardas; (f) Tavares; and (g) South. The gray area corresponds to the 95% highest posterior density limits for the effective population size. The thin and thick dotted lines are the lowest and the mean estimated tMRCA. Below the plot is the time represented in years, and in the left side are the effective population sizes.

The Mantel’s test detected a positive and significant correlation between genetic and geographical distances, for the overall geographical range in CC mtDNA (r=0.38, P=0.00), microsatellite data (r=0.53, P=0.00) and for chromosomal rearrangements (r=0.60, P=0.00). The partial Mantel tests do not detect significant correlations between genetic and chromosomal variations, while partial out the geographical distance, in both CC mtDNA (r=−0.07, P=0.73) and microsatellite data sets (r=0.10, P=0.20).

Clear positive phylogeographic structures were also observed through the spatial autocorrelation analyses (Supplementary Figure S5). The mean genetic distances (Ay) estimated for CC mtDNA, microsatellites and chromosomal rearrangements were 0.015, 0.770 and 0.363, respectively. For CC mtDNA data, highly positive correlations were evident among the first three distance classes (from 0 to 42 km). For microsatellite data, the positive structure persisted up to around 56 km, and for chromosomal rearrangements, the strongest correlations were observed up to around 70 km.

Population structure

The microsatellite pairwise estimates of F_ST and R_ST revealed moderate to high levels of genetic differentiation among sampling sites, ranging from 0.09 to 0.61 for F_ST, and from 0.05 to 0.79 for R_ST. Most pairwise comparisons showed values higher than 0.3 (Supplementary Table S4).

The AMOVA results for the three mtDNA data sets, using eight groups of parental karyotypes showed most of the genetic variation distributed among populations within karyotypes (49.75–62.16%). On the other hand, the test in which the populations were subdivided into seven groups, considering present and past geographical barriers, revealed higher percentages of genetic variation at the group level, ranging from 55.58 to 58.28% (Table 3). For microsatellite data, the two AMOVA tests showed more than 67% of the genetic variation within populations; however, the percentage of variation among karyotypic groups was slightly higher when compared with geographical groups, 10.49% versus 8.77%, respectively (Table 3).

Table 3 Results of analyses of molecular variance for mitochondrial DNA data sets and microsatellites

The most likely number of clusters in the Structure analysis, according to the ΔK approach of Evanno was 12 (Supplementary Figure S6). The major divergences among the 60 runs at K=12 were regarding the estimated membership coefficients (q) and the assignment of specimens from the localities of MC and SJN. The samples from MC showed a private cluster in some runs, whereas in other runs, the specimens showed mixed ancestries (q<0.8) between the clusters representing the neighboring localities of JG/ILH and GAI. The locality of SJN kept all the individuals in a private cluster in some runs, and in other runs, the specimens were grouped in the same cluster of the neighboring localities of BJ1 and BJ2. All the other sampling sites showed consistent results among different runs. The plot of population clustering retrieved in the run with the highest ln Pr (X/K)=−13 404.3 is represented in Figure 5a. In the Tess analysis, the plot of DIC × K_max began to approach a plateau around K=12 (Supplementary Figure S7). We analyzed the best runs for Ks from 12 to 16. All these Ks retrieved the same pattern of grouping localities in 12 main clusters, and the cluster retrieved above 12 comprised only minor subdivisions in the sampling sites LLB, FOR, 53, 64, TV1 or TV2 (data not shown).

Figure 5

Bayesian clustering and specimen assignments for the 12 clusters identified by the (a) Structure and (b) Tess programs. Each specimen is represented by a single bar, and each cluster by a color. Sampling sites are plotted geographically in the North–South direction, from left to right. Above plot: parental parapatric karyotypes, hybrid zones among them (HZ) and main possible geographical barriers. Below plot: BAR, Barco Beach; BJ1, Bujuru 1; BJ2, Bujuru 2; EBL, East Barros Lake; FOR, Fortaleza Lake; FSM, Farol de Santa Marta; GAI, Gaivota Beach; GUA, Guarita Beach; ILH, Ilhas; JG, Jaguaruna; MC, Morro dos Conventos; MOS, Mostardas; PAL, Palmares do Sul; PAS, Passinhos; PIT, Pitangueira; PST, Passo de Torres; OSO, Osório; 35, road km 35; 53, road km 53; 64, road km 64; 96, road km 96; 108, road km 108; 115, road km 115; SBL, South Barros Lake; 17S, 17 km south of Mostardas; 26S, 26 km south of Mostardas; SJN, São José do Norte; TRA, Tramandaí; TV1, Tavares 1; TV2, Tavares 2. *Roman numbers correspond to Structure clusters, and the uppercase letters correspond to Tess clusters.

The 12 clusters retrieved by both Structure and Tess tended to group together specimens from the same sampling site and nearby localities, in a similar pattern. One of the main discrepancies was with respect to MC, which was completely admixed for Structure among several clusters, but in Tess its ancestry came from a single cluster (B; Figure 5b). Moreover, Tess showed a smaller proportion of individuals with admixed ancestries (q⩽0.8) compared with Structure (∼10% versus 27%, respectively). This difference is mainly associated with localities MC, BAR, TRA, OSO, PAS, PIT, 35, 53, 64, 96 and 115, in which most specimens could be classified as non-admixed by Tess, but had mixed ancestries in the Structure results. On the other hand, sampling site 108 showed more admixed specimens in Tess than in Structure.

Discussion

Isolation by distance

C. minutus has a nearly linear geographical distribution along the southern Brazilian coastal plain. In general, when species of small mammals have a linear distribution, dispersal and gene flow will be geographically limited to neighboring populations, mainly because of their low ability to disperse over significant distances relative to the species’ range, which is a common feature in subterranean rodents (Busch et al., 2000) and may result in a pattern of isolation by distance. Phylogeographical studies based on control region, COI and/or cytochrome b mtDNA data revealed that species of tuco-tucos such as C. pearsoni, C. flamarioni, C. australis and C. lami have low levels of gene flow among populations, and a pattern of isolation by distance associated with their narrow geographical distribution (Mora et al., 2006; Fernández-Stolz, 2007; Tomasco and Lessa, 2007; Lopes and Freitas, 2012).

C. minutus showed a geographical limitation on the distribution of mtDNA haplotypes, as most of them were limited to single sampling sites, and shared haplotypes were mainly among neighboring localities (Figures 1 and 3a and Table 1). The most widely separated localities that shared haplotypes, for CR, COI and CC data, respectively, were FOR and 35 (41.5 km apart, haplotype CR20), SBL and MOS (140 km apart, haplotype COI8) and 96 and MOS (32 km apart, CC13; Supplementary Table S3). Bayesian clustering of populations based on microsatellite data tended to group together specimens from the same sampling site and nearby localities (Figures 5a and b). Likewise, the distribution of parental karyotypes was parapatric, and no specimen was found outside the limits of its source karyotypic population (Table 1).

Analyses of the variation in the skull shape of C. minutus found a weak but significant association between morphometric and geographical distances, and the authors suggested that differences among populations follow a subtle pattern of isolation by distance (Fornel et al., 2010).

The Mantel tests and the spatial autocorrelation analyses showed overall positive and significant correlations between genetic and geographical distances, for the mitochondrial, nuclear and chromosomal data sets, and the strongest positive structures were observed up to around 70 km. However, the values of Ay did not increase linearly with increasing geographical distance, showing an unstable pattern of variation for all data sets, mainly from the 10th distance class, corresponding to 140 km onwards (Supplementary Figures S5A–C). This pattern was described by Diniz-Filho and Telles (2002) as a ‘stabilizing profile’, which may reflect influences of positive genetic/geographical associations among groups of closely distributed populations, and phylogeographical groups of populations primarily structured by other historical events, for example, the genetic structures revealed by the seven haplogroups for mtDNA data or the 12 clusters of localities retrieved in the Bayesian clustering analyses. The sharp increase in genetic and geographic distance correlations, until it reaches an irregular pattern, is expected under a model of isolation with asymmetric gene flow in a metapopulation. A similar pattern was observed in the perrensi group of ctenomyids (Fernández et al., 2012).

Geographical barriers

The phylogenetic tree generated from concatenated data (Figure 3a) showed a basal polytomy between the C. minutus branches and the outgroup. Polytomies are commonly found among and within species of the genus Ctenomys, and have been suggested to be the result of an early simultaneous cladogenesis event (Lessa and Cook, 1998). Nevertheless, we were able to recover a pattern of genetic structure of the sampling sites subdivided in seven main haplogroups.

Time estimates and calibration of molecular clock is often a hard task because of the scarcity and uncertainty of the fossil records for some taxa, sometimes resulting in substitution rates with large confidence intervals, which may account for the largely overlapping 95% HPD of tMRCA calculated for the seven mtDNA haplogroups of C. minutus. Nevertheless, we can observe a general trend of older mean estimates of tMRCA for the haplogroups of C. minutus came from its northern range, with exception of North 1, and toward the south the mean tMRCA decreased, reaching its minimum in the South haplogroup (Figure 4). Moreover, the southern locality of SJN has the lowest levels of microsatellite variability, showing a total of six monomorphic loci. This pattern may result from a more ancient genetic structure achieved in the northern geographical distribution, while the southern sampling sites may account to a founder effect of more recent occurrence.

Around 200 ka, the sea level started a regressive cycle, and in 141 ka, the southern Brazilian coastal plain became wider than nowadays. During this period, both Paleojacuí north and Paleocamaquã were present in the region. The subsequent transgressive cycle, around 133 ka, closed completely the Palejacuí north outlet to the sea. A new glacial cycle led to a long and oscillating regressive phase of the Atlantic sea level, which persisted for around 100 000 years and reached its maximum at 18 ka. The Paleojacuí south arose on the coastal plain in the beginning of this regressive cycle. Both the Paleojacuí south and the Paleocamaquã were completely closed by sediments 6 ka, during the last great transgressive cycle of the Atlantic sea level (Tomazelli et al., 2000; Weschenfelder et al., 2008; Baitelli, 2012).

These three paleochannels may have reduced the gene flow between populations from opposite sides of the rivers during the Pleistocene. As a consequence, we detected associated patterns of genetic structure that have persisted to the present, which can be observed in Bayesian clustering analyses (Figure 5) and in mtDNA haplogroups (Figures 1 and 3a and b). The persistence of these genetic structures may have been reinforced by the metapopulation structure, common to ctenomyids, in which the effects of low levels of gene flow hamper the gene exchange between the subpopulations, and the genetic drift associated with the small subpopulation sizes maintain the low levels of genetic diversity.

The transition between the sand fields and dune line represent another environmental discontinuity along the range of C. minutus, which may inhibit free dispersal among specimens from different habitats. For mtDNA data, we recovered the Coastal haplogroup (corresponding to specimens from the dune line) and Barros Lake (sand fields; Figure 1). For microsatellite data, a pattern of genetic structuring following the transition between these two environments was less pronounced. The tuco-tucos living near the transition region showed admixed ancestries (q⩽0.8) among one main cluster (IV–D for specimens from the dune line, and V–E for specimens from sand fields), and several other less-representative clusters (Figure 5).

Environmental characteristics such as soil hardness lead to morphological adaptations in tuco-tucos, mainly in the anterior and posterior limbs and in the cranium, as they are burrowing rodents and use their teeth, paws and legs to dig tunnels. Fornel et al. (2010) found significant differences in skull morphology between specimens of C. minutus that live in sand fields or dunes. Therefore, these different environments may also be associated with the genetic discontinuities observed in our data.

The role of the Mampituba and Araranguá Rivers as effective geographical barriers to gene flow were conflicting regarding mitochondrial and nuclear molecular markers. One possible explanation for the absence of clear genetic structuring caused by the Mampituba and Araranguá Rivers may lie in the past transient status of the river mouths (Tomazelli and Villwock, 2000; S Dillenburg, personal communication), which may have often allowed gene flow between populations from opposite banks, which have not reached reproductive isolation. Nowadays, the banks of both rivers have stabilized, hampering the contact between specimens from opposite sides.

Most of mtDNA genetic variation is apportioned among groups geographically delimited, according to AMOVA (Table 3), but microsatellite data showed the lowest proportions of genetic variation among groups, suggesting the involvement of factors other than geographical barriers or chromosomal rearrangements affecting the partition of genetic variation in nuclear genome. The potential homoplasy and high levels of polymorphism in microsatellites, and the lower effective population sizes of mtDNA, may account for different results observed in AMOVA tests.

Cytonuclear genome discordance has been described for some vertebrates at both intra- and interspecific levels. Genome components, with distinct modes of inheritance, commonly show differences in phylogeographic patterns and levels of introgression between hybridizing species/populations. These differences can result from distinct mutation rates and effective population sizes among molecular markers, and also from sex-biased dispersal (Petit and Excoffier, 2009). The mtDNA has non-recombinant maternal inheritance, which allows the retention of haplotypes shared between populations recently isolated, or to keep the differentiation between populations in recent secondary contact that had previously differentiated in allopatry. Whereas, male-biased dispersal, that is characteristic of the species of the genus (Malizia and Busch, 1991; Cutrera et al., 2006), can lead to faster and farther dissemination of alleles of nuclear molecular markers relative to mtDNA, and associated with the nuclear genome recombination and high mutation rates of microsatellite loci, it can mask the signals of historical evolutionary processes in the nuclear genome (Petit and Excoffier, 2009; Yang and Kenagy, 2009).

The skull morphometric data for C. minutus indicate that there was no direct relationship between Robertsonian fusions/fissions and skull morphology, and these characters are independent, although both vary with the geographical distribution (Fornel et al., 2010). A phenotypic structuring into three morphological groups was suggested for C. minutus: (i) the populations in the northern part of the geographical range, isolated from the other localities by the Araranguá River; (ii) populations in the southern part of the distribution, isolated by the Paleojacuí south, which might be reinforced by the pericentric inversion; (iii) and all other populations associated with lakes and lagoons, in the middle of the distribution.

The geographical features of the southern Brazilian coastal plain may also be associated with the fixation of new chromosomal rearrangements in C. minutus, and two different scenarios of a step-by-step mechanism of chromosomal evolution can be suggested for this species. The first hypothesis relies on the metapopulation structure, with small subpopulations distributed in patches of suitable habitats (Fernández et al., 2012). Whether a new chromosomal rearrangement causes little or no negative meiotic effects, it can be more easily fixed in peripheral and/or small isolated populations (Baker and Bickham, 1986) by means of genetic drift, spreading to new, recently colonized areas. The second hypothesis is an extension of the first, also based on the fixation of chromosomal rearrangements by means of genetic drift, but in this case the chromosomal rearrangements remain isolated from each other in allopatry because of effective geographical barriers to gene flow. An example from this study is the transition between karyotypes 2n=48c and 46a, which are delimited by the Araranguá River (Figure 1). As these barriers no longer exist, the gene flow between specimens from parapatric karyotypes, which were not reproductively isolated, could give rise to secondary-contact hybrid zones, for example, between karyotypes 2n=42 × 46b, overlapping the region corresponding to Paleojacuí south (Figure 1).

Chromosomal rearrangements

Subterranean rodents show some of the highest rates of chromosomal evolution known in mammals (Reig et al., 1990). The cytogenetic data available for C. minutus reveal a total of 45 distinct karyotypes (Freitas, 1997; Gava and Freitas, 2002, 2003; Freygang et al., 2004; Castilho et al., 2012; this study), and most of this polymorphism was retrieved in the intraspecific karyotypic hybrid zones.

The six hybrid zones identified for C. minutus arose even between parapatric karyotypes differing by one Robertsonian rearrangement: (i) 2n=50a × 48c; (ii) 46a × 48a; (iii) 46b × 48b; (iv) and 48b × 50b; or more than one: (v) 2n=48a × 42; and between karyotypes differing by Robertsonian rearrangements associated with a pericentric inversion: (vi) 2n=42 × 46b (Figure 2). Adaptability and fertility in chromosomally heterozygous carriers can be maintained by mechanisms that partially or totally suppress recombinations with harmful effects (White, 1978; Rieseberg, 2001). Heterozygotes for simple Robertsonian centric fusion/fission rearrangements, for example, can produce balanced meiotic systems, by means of trivalents, which often show normal segregation during meiosis (Baker and Bickham, 1986).

Different patterns of genetic structures and levels of introgression were observed among the hybrid zones of C. minutus. These patterns may vary depending on the molecular marker analyzed, which include differences in mutation rates, recombination and effective population sizes; the kind of chromosomal rearrangements involved and whether they are implicated or not in reductions in gene flow between chromosomally divergent populations; whether there are different patterns of crossings, or both sexes from both parental karyotypes are able to mate and breed equally; or whether the hybrid zone is associated or not with a geographical barrier and the effectiveness and persistence of this geographical feature as a barrier to gene flow between populations.

The recent literature about ctenomyids provides some examples of karyotypic polymorphic species that fail to show genetic structure following a pattern of chromosomal rearrangements (Lopes and Freitas, 2012). For C. pearsoni, it is suggested that new chromosomal rearrangements, besides being a common feature, seem not to be followed by sterility, reductions in fitness or negative heterosis of heterozygote carriers (Tomasco and Lessa, 2007). Studies by Gava and Freitas (2002, 2003), using cytogenetic data for the C. minutus hybrid zones between 2n=46a × 48a and 2n=48a × 42 showed no evidence of underdominance, and demonstrated that the chromosomal polymorphisms fail to cause sterility in heterotypes. The chromosomal rearrangements in the ‘perrensi’ group were considered to be the result of an ongoing isolation with a migration process, rather than representing a barrier to gene flow (Caraballo et al., 2012; Fernández et al., 2012). Some highly chromosomally polymorphic mammals, such as Mus musculus domesticus and Sorex araneus, show a step-by-step mechanism of karyotype evolution, similar to that of C. minutus. However, although the presence of karyotypic populations is well defined, molecular markers do not match this subdivision, and the karyotypic populations are still considered as part of a single species (Wójcik et al., 2002; Piálek et al., 2005).

C. minutus shows no chromosomal polymorphism within the parental karyotypic populations (outside hybrid zones), and differences among neighboring parental karyotypes, although minor, are fixed. However, chromosomal rearrangements are not strongly associated with the patterns of nuclear genetic structure observed, as suggested by AMOVA, partial Mantel test and Bayesian clustering analyses, and thus probably has no or very little role in the diversification of C. minutus populations. As suggested by Lukhtanova et al. (2011) ‘we cannot exclude that geographically distant and chromosomally divergent populations would display reduced fertility if crossed, although they are connected by a chain of compatible populations that should allow gene flow’.

In summary, the effects of genetic drift, low levels of gene flow relative to the geographical range, the metapopulation structure and the geographical features of the southern Brazilian coastal plain may account for the patterns of genetic/geographical structure and the step-by-step mechanism of chromosomal evolution observed in this study for C. minutus. Further researches regarding the chromosomal rearrangements and their effects on heterozygote carriers will allow to better understand the evolutionary process in which C. minutus is involved.

Data archiving

Sequence data have been submitted to GenBank. Ctenomys minutus accession numbers: HM236969–HM237008 and HM237009–HM237043; Ctenomys torquatus: HM443438 and HM443439; Ctenomys flamarioni: JQ341041 and JQ341052. Microsatellite data have been submitted to Dryad repository, doi:10.5061/dryad.d12j8.