A comprehensive overview of the genetic diversity in Thylamys elegans (Didelphimorphia: Didelphidae): establishing the phylogeographic determinants

Abstract

Background

For the genus Thylamys, the rivers have been reported as barriers to dispersal, limiting current and historical distribution of its lineages. We hypothesized that the Maipo river has affected the genetic structure of northern and southern lineages of Thylamys elegans, recovering a phylogenetic relationships with reciprocally monophyletic sister groups on opposite river banks. We evaluated the role of other rivers in the Mediterranean zone of Chile as historical and recent modulators of the biogeographic processes of this species.

Methods

We applied a phylogeographic approach, using the cytochrome-b mitochondrial gene for 93 individuals of T. elegans, from 37 localities in a latitudinal gradient between 21°25’ and 35˚56’S, encompassing a geographic area between the Atacama Desert and most of the Mediterranean Chilean zone.

Results

The phylogenetics results recovered six lineages within T. elegans: Thylamys elegans elegans, Thylamys elegans coquimbensis, the Loa lineage and three other lineages not described previously (Aconcagua, South 1 and South 2). We suggest that following rivers play a role like primary barrier: the Maipo river in the genetic differentiation of northern and southern ancestral lineages, and the Mataquito river and its tributary Teno river for the South 1 and South 2 lineages. On the other hand, the Quilimarí river preserve the genetic divergence in T. e. coquimbensis and Aconcagua lineage and the Aconcagua river in Aconcagua lineage and T. e. elegans acting like secondary barriers.

Conclusions

We concluded that the genetic diversity and biogeographic history of T. elegans was shaped by mountain glaciers, changes in river water levels during the Pleistocene glaciations and hyperaridity, promoting the differentiation and persistance of the T. elegans lineages.

Background

Since the 19th century, rivers have been proposed as geographical barriers to dispersal for several taxa [1]. The riverine barrier hypothesis proposes that populations of different river banks become progressively isolated each other and predicts populations of opposite banks with reciprocally monophyletic relationships [2, 3]. The effects of rivers on the biota has become more evident by the use of molecular and bioinformatics tools. Three hypotheses are proposed to explain riverine divergence by mean of primary differentiation, secondary contact and dispersal: (1) primary differentiation, with a river acting as a complete vicariant barrier for an existing taxon range, resulting in a topology with reciprocally monophyly of sister clades from opposite banks; (2) secondary differentiation, with a river as a common contact zone acting as a barrier to dispersal, preventing secondary contact of previously differentiated populations and a topology with opposite river lineages showing more recent common ancestor with others populations of the same river bank and lineages from opposite banks are not sister clades; (3) dispersal events from an established populations to the opposite riverbank, with the river being a permeable barrier to gene flow and a topology with paraphyly [2,3,4]. If rivers represent a barrier to populations of some species, it is expected that the diversification between populations from opposite riverbanks coincide with the formation of the river. If they do not coincide temporarily, the river is unlikely to be the primary divergence agent, but then it could have acted as a barrier to gene flow [5].

Didelphidae is one of the oldest families of extant mammal [6]. Species of this family inhabit in forest and open habitat in the Neotropics [7]. Rivers have been reported as geographic range boundaries or as barriers to gene flow in some populations of didelphid marsupials of genus Philander, Marmosa, Gracilinanus, Marmosops and Cryptonanus in the Amazon River [8,9,10]; the Paraná and Paraguay rivers [11, 12], the São Francisco river [13, 14] the Paraguay river [15], the Tocantins river [14], the Doce and Paraíba do Sol rivers [16] and Tapajós rivers [10]. In Thylamys, the effect of rivers such as geographic barrier to gene flow or geographic boundary of distributional range has been reported for the following rivers: Paraná for T. citellus and T. pulchellus; Bermejo for T. pulchellus and T. pusillus, Paraguay for T. pusillus, Paraná and Uruguay for T. citellus [17], Pilcomayo for T. pusillus A and T. pusillus B; and Paraná and Uruguay for T. pusillus C [18]. Paraná for T. pusillus C and T. pusillus A, and Paraguay for T. pusillus B, T. pusillus C and T. pusillus A [18], and Cañete for T. tatei and T. sp [19].

Thylamys elegans is one out of two Thylamys species present in Chile [19, 20]. This species has a huge distribution from the north of the Loa river Mouth (21°25’S) to Angol and adjacent areas of the Nahuelbuta coastal range (37°50’S) according to historical records, covering an extension close to 1850 km [19,20,21,22]. Thylamys elegans  is a mouse-like marsupial whose activity is mainly arboreal and crepuscular. Its diet items are composed mainly by insects although fruits, small vertebrates and even carrion occasionally are consumed [7, 23]. This and other species of Thylamys genus are characterized by its tail incrassation during winter periods when the specimens fall in torpor [23,24,25]. A marked north-south climatic gradient is present in this wide distribution, with different environment in each geographical zone [26,27,28]. Thus, those populations that inhabit xeric environments are associated with small patches of available vegetation (see Fig. 1b), while the specimens that inhabit in humid environment had several different habitats available even into brushy ecotones of relictual cloud forest of the Coquimbo region, northern and central part of the Coastal ranges and the steppes of Central Chile [19, 20, 23, 29, 30] (See Fig. 1c). Clearly this huge and contrast distribution produce a strong cline in size, color and shape of the different populations [20, 21, 23, 24, 31, 32]. Actually, the differences between clines have been enough to suggest three subspecies based only in morphological characters [21, 23, 24, 32], but its distribution boundaries are still uncertain [19, 20, 33].

Fig. 1

Morphological variability and environment of Thylamys elegans. a Adult specimen and landscape from the Loa river Mouth, Loa lineage (Clade A); b Adult specimen corresponding to T. e. coquimbensis from Paihuano and ravine from Llanos de Challe (Clade B); c Adult specimen and landscape from Los Molles corresponding to the Aconcagua lineage (Clade C); d Adult specimen corresponding to T. e. elegans from Punta Curaumilla and landscape from Las Docas (Clade D; e) Adult specimen from Lipimávida and landscape from Pupuya, corresponding to the South 1 lineage (Clade E); f Adult specimen (not sequenced) and landscape from Nahuelbuta Coastal range, probably corresponding to the South 2 lineage (Clade F). Northern clade encompasses the Loa, T. e. coquimbensis, Aconcagua and T. e. elegans lineages and southern clade encompasses the South 1 and South 2 lineages

Additionally, studies have found genetic structure in the distribution of T. elegans [19, 20, 34], but the genetic clades not match with the morphological subspecies described (see [20, 33] for T. e. coquimbensis). These evidences suggest that morphological variability could be more related with the environment than the genetic fingerprint and for that reason is still difficult delimit the distributional range of the different lineages or subspecies. Anyway, so far, what is limiting the distribution of each lineage or subspecies in T. elegans is still uncertain. A previous study proposed in Mediterranean zone of Chile, the Maipo river as an strong geographic barrier in the genetic differentiation of northern and southern clades of Thylamys elegans [19] and a genetic differentiation congruent with the position of the Maipo river has been reported for other taxa [35,36,37,38,39,40,41]. Indeed, some rivers can act as barriers to gene flow and/or geographic boundary of distribution structuring genetic diversity in different taxa from Chile [36, 41,42,43,44,45,46]. However, none of these papers discussed in detail if the riverine divergence would be a primary or secondary differentiation [47]. On the other hand, the historical and recent effects of other rivers different to Maipo river on the differentiation of T. elegans lineages are still unknown, being necessary to evaluate if the role of rivers was for primary diversification, secondary contact or dispersal [4]. Knowing these effects is very relevant in the geographical context in which this species is distributed, where the spatial complexity of the landscape is proposed as the cause of the high endemism of one of South America’s biodiversity hotspots [48, 49]. Therefore, under these evidences the hypothesis about the rivers as phylogeographic determinant is novel to test it in T. elegans given that in Chile a role of rivers in the genetic differentiation of a mammal has never been discussed in detail.

Although the largest rivers in South America are in tropical and subtropical regions, the effects of smaller rivers in the temperate region as geographic boundaries or barriers to gene flow were reported in Chile for the Huasco, Copiapó, Elqui, Aconcagua, Maipo, Mataquito, Claro, Lircay, Maule, Itata and Bío-Bío rivers [35,36,37,38,39,40,41,42,43,44,45,46, 50,51,52,53,54,55]. Therefore, the main goal of this study was to evaluate, through a phylogeographic approach, the role of some rivers of Mediterranean Chilean zone on the genetic structure of T. elegans. We propose that some riverine barriers from Mediterranean zone of Chile shape the genetic diversity of Thylamys elegans, triggering the genetic divergence of its lineages through vicariant events as primary riverine divergence, or persistence of its genetic differentiation with rivers acting as geographic boundaries of lineages previously differentiated and posterior secondary riverine differentiation.

Methods

Collections

We perform nine field campaigns between 2010 and 2012 which covered most of the distributional range of T. elegans. From these campaigns and voucher specimens stored at the Colección de Flora y Fauna Patricio Sánchez Reyes, Pontificia Universidad Católica de Chile, we had available samples from north of the Loa river Mouth, (21°25’S) in the Tarapacá region to Tregualemu (35°56’S) in the Maule region (Fig. 2, see details in Additional file 1: Appendix I). Specimens were collected with Sherman traps (8 × 9 × 23 cm) using a mixture of oat, banana and canned fish as bait or crushed oat. Captured individuals were euthanized with isofluorane overdose and cervical dislocation. Collection and handling recommendations were followed according to the American Society of Mammalogists [56]. Samples were obtained from different tissues such as blood, hair or liver. Also from ear skin or tail from corpses found on fieldwork and donations of tissues by collaborators. Eighteen samples correspond to donations, 63 samples belonged to Colección de Flora y Fauna Patricio Sánchez Reyes, Pontificia Universidad Católica de Chile and 14 samples from 10 localities were captured and euthanized in this study (See details in Additional file 1: Appendix I). The samples were stored in cryotubes preserved in liquid nitrogen or 70 % ethanol. All voucher specimens (tissues, skins or skeletons) were stored at the Colección de Flora y Fauna Patricio Sánchez Reyes, Pontificia Universidad Católica de Chile.

Fig. 2

Summarized phylogeny and map with the 37 sampled spots in this study. The colors of the sampled spots match with clades on the tree (also for Fig. S1 and Fig. S2). Most important branches of the rivers that limit the lineages of T. elegans are shown (Quilimarí, Aconcagua, Maipo, Mataquito and its tributary Teno). Loa lineage (Clade A: blue); T. e. coquimbensis (Clade B: cyan); Aconcagua lineage (Clade C: orange); T. e. elegans (Clade D: brown); South 1 lineage (Clade E: green); South 2 lineage (Clade F: red). Northern clade encompasses the Loa, T. e. coquimbensis, Aconcagua, and T. e. elegans lineages and southern clade encompasses the South 1 and South 2 lineages

DNA extraction, amplification and sequencing

DNA was extracted with phenol-chloroform procedures [57]. The complete cyt-b mitochondrial gene (1149 base pairs) was used for 93 individuals of T. elegans from 37 localities and 2 sequences corresponding to T. pallidior A and T. pallidior B respectively as outgroup, (including 32 sequences from previous studies in Thylamys genus [19, 20, 58] and 63 new sequences generated in this study with the following genbank access number (MZ868645- MZ868707, see details in Additional file 1: Appendix I). This gene has been widely used in phylogeographic studies in Thylamys genus because its variability [17, 20, 33, 58]. PCR amplification reactions followed described protocols [19, 20], and products were purified with QIAquick PCR Purification Kit (QUIAGENTM Inc., Valencia, CA, USA), and sequenced in MACROGEN Inc. (Seoul, Korea). Sequences were edited in Bioedit [59] and aligned in Clustal W [60] using default parameters and encoded to aminoacid in DnaSP v6 [61] to verify the lack of internal stop codons.

Phylogenetic reconstruction

The phylogenetic reconstructions was based on 49 haplotypes and 95 sequences using Maximum Likelihood (ML) criterion on the online version of IQ-TREE software [62] available at http://iqtree.cibiv.univie.ac.at [63]. The best model of evolution was retrieved with ModelFinder [64] implemented in IQ-TREE using the Bayesian Information Criterion (BIC) prior to the construction of the ML trees. To assess branch support, the ultrafast bootstrap approximation (UFboot2) with 1000 replicates was implemented [65]. Bayesian Inference (IB) was performed using Mr. Bayes 3.2.6 [66]. Four simultaneous Markov chain Monte Carlo chains were run for 10,000,000 generations, with a sampling frequency of 1,000 generations. The first 25 % trees were discarded as burn-in and remaining trees were used to compute a 50 % majority rule consensus tree and to obtain posterior probability estimates for each clade. The ML and BI phylogenies were rooted using T. pallidior A and T. pallidior B as outgroup (see details in Additional file 1: Appendix I, Fig S1, Fig S2) following to [19, 58]. To assess the genetic distance between lineages within T. elegans a corrected distance Kimura 2-parameters (K2P) was calculated in MEGA 7.0 [67].

Genetic structure

To estimate the geographic population structure of T. elegans in a Bayesian framework, we used GENELAND v4.0 [68]. The number of genetic clusters (K) was determined using a prior of K between 1 and 37 (number of localities used in this study). We performed 10 independent runs of 5,000,000 iterations, sampling every 1,000 steps. Burning and convergence of the chains were determined with a 0.1 cutoff value and using uncorrelated and correlated frequency models, selecting the most probable model by using the Bayes Factor statistic [69]. Correlated frequency model was chosen and we set the K value to 6 considering the posterior probability distribution of the parameter. After that, five runs were post-processed by burning 5 % of first iterations and we obtain the posterior probability of assignment for each individual to a genetic cluster and to each pixel in the spatial domain. Maximum number of Poisson-Voronoi tessellation nuclei was set at 279 according to the number of genotypes used for the analysis (93 × 3) [70].

Divergence times

To estimate the divergence times for the T. elegans lineages, we used the Bayesian algorithm implemented in BEAST 2.5.2 [71]. The best molecular clock model was selected by Bayes Factor comparing three clock models (strict, uncorrelated lognormal relaxed and uncorrelated exponential relaxed) in Tracer v1.5 [72]. The uncorrelated lognormal relaxed molecular was chosen, using as priors the general time reversible model of sequence evolution (GTR + G) and a Yule speciation process. The calibration points selected for the analysis following to [19]were the same T. pallidior-T. elegans molecular divergence, estimated at 6.11 Million year ago (Mya), (Highest posterior density HPD: 5.92-6.3), and the T. tatei-T. sp molecular divergence at 4.65 Mya (HPD: 4.46-4.85) (Fig S3). Estimate of posterior distribution was obtained through Markov Chain Monte Carlo MCMC method with 20,000,000 iterations, sampling parameters every 10,000th steps, and burning the first 10,000 to achieve convergence of posterior probability distribution. The convergence of samples obtained by MCMC method was visualized in Tracer v1.5 [72]. We verified if the Effective Sample Size of all parameters had values greater than 200. Maximum credibility tree was calculated in TreeAnotator v1.7.0 [73]. The tree was displayed and edited in Figtree 1.4.2 [74].

Results

BIC identified the HKY+F+G4 model as the best nucleotide substitution model for a matrix with 49 haplotypes and other with 95 sequences. Gamma shape distribution parameter was 0.2818 for haplotype matrix and 0.183 for sequences matrix. Both ML and IB reconstruction with haplotypes and all sequences showed the same phylogenetic relationship among six lineages (see Fig. 2, Fig. S1 and Fig. S2). T. elegans is a clade highly supported (100/1). The ingroup exhibited six monophyletic groups with high bootstrap supports and posterior probability values (Fig. S1 and Fig. S2). The southern clade contained the South 1 and South 2 lineages (99/1), and the northen clade contained T. e. elegans sister to (T. e. coquimbensis, (Loa lineage, Aconcagua lineage)) with high support (97/1). The Loa lineage (Clade A) is restricted to north of the Loa river mouth (21°25’S) (1 locality, n = 6). The Aconcagua lineage (Clade C) encompasses localities from Fundo el Roble (32°16’S) to Fundo Chuico Blanco (32°48’S) covered about 100 km from the Quilimarí river and the Santa Inés-Imán mountain range to the Aconcagua river and the Chacabuco range (7 localities, n = 15). T. e. coquimbensis (Clade B) encompasses localities from Quebrada El León (26°57’S) to Valle El Mauro (31°58’S). The southern boundary would be the Quilimarí river and Santa Ines-Imán mountain range (8 localities, n = 19). On the other hand, T. e. elegans (Clade D), encompasses localities from La Campana (32°57’S) to Viña Leyda (33°34’S) between the Aconcagua river and Chacabuco range and the Maipo river and Altos de Cantillana range (10 localities, n = 33). The South 1 lineage (Clade E) encompasses localities from San Enrique (33°53’S) to Duao (34°52’S), between the Maipo to Mataquito river and its tributary Teno river (7 localities, n = 13). The South 2 lineage (Clade F) encompasses localities from Los Queñes (35°0’S) to Tregualemu (35°56’S), between the Mataquito and Teno rivers to Tregualemu (4 localities, n = 7) (Fig. S1, Fig. S2).

The corrected distance (K2P) values were between 2 % for the South 1 and South 2 lineages to 11.9 % for the Loa and South 1 lineages (Table 1). On the other hand, our results show a genetic distance of 5.5 % between the northern clade (i.e. Loa lineage, T. e. coquimbensis. Aconcagua lineage and T. e. elegans) and southern clade (i.e. South 1 lineage and South 2 lineage).

Table 1 Corrected genetic distance values between the six lineages of T. elegans using the K2P model

GENELAND analysis detected six genetic clusters (Fig. 3) which fit with the six lineages of T. elegans recovered in the phylogenetics reconstruction (Fig. 2, Fig. S1 and Fig. S2). Detection of the geographical genetic discontinuities correspond to the Loa, Quilimarí, Aconcagua, Maipo, Mataquito and Teno rivers and the Santa Inés-Imán, Chacabuco and Altos de Cantillana mountain ranges previously described (Figs. 2 and 3). GENELAND (Fig. 3) estimated clusters corresponding to the following lineages described in Fig. 2, Fig. S1 and Fig. S2: (1) Loa (Clade A); (2) T. e. coquimbensis (Clade B); (3) Aconcagua (Clade C); (4) T. e. elegans (Clade D); (5) South 1 (Clade E), and (6) South 2 (Clade F). T. e. coquimbensis encompass 550 km, and the Aconcagua, T. e. elegans, South 1 and South 2 lineages together encompass 365 km (Fig. 3).

Fig. 3

Genetic structure in geographic space of T. elegans. Maps of posterior probabilities isoclines of genetic clusters membership with its respective populations belonging to a panmictic unit correspond to the six lineages of T. elegans detected by GENELAND analyses. Last column: synthetic map of the K-origin of each pixel. Blue line is the Maipo river such a geographic boundary between the northern and southern clades; other six columns: relative posterior probability of belonging to each K cluster. Black dots show sampled populations and lighter background colours indicate higher posterior probability of membership to the given K cluster. Thinner black lines are isoclines of inclusion probability and numbers are posterior probabilities values. Cluster 1 correspond to Loa lineage (Clade A), Cluster 2 correspond to T. e. coquimbensis (Clade B), Cluster 3 correspond to Aconcagua lineage (Clade C), Cluster 4 correspond to T. e. elegans (Clade D), Cluster 5 correspond to South 1 lineage (Clade E), and Cluster 6 correspond to South 2 lineage (Clade F). Northern clade encompasses the Loa, T. e. coquimbensis, Aconcagua and T. e. elegans lineages and southern clade encompasses the South 1 and South 2 lineages

The T. elegans clades diverged as follow (see Fig. 4 and Fig. S3): northern and southern clades, 3.14 Mya (HPD: 2.95-3.33); T. e. elegans from the clade with T. e. coquimbensis sister to the Loa and Aconcagua lineages, 1.84 Mya (HPD: 1.41-2.27); T. e. coquimbensis from the Loa and Aconcagua lineages, 1.41 Mya (HPD: 1.07-1.8); the Loa and Aconcagua lineages, 1.28 Mya (HPD: 0.91-1.64); the South 1 and South 2 lineages, 0.81 Mya (HPD: 0.48-1.2) (Fig. 4); the Andean-Coastal clade from T. pallidior A and T. pallidior B, 3.35 Mya (HPD: 3.16-3.54), T. sp and T. tatei, 4.65 Mya (HPD: 4.46-4.85) and T. pallidior-T. elegans from T. tatei-T. sp, 7.64 Mya (HPD: 6.61-8.73) (Fig. S3).

Fig. 4

Divergence times of the six lineages of T. elegans. Clades correspond to Fig. 2. Number of nodes are average values and bars in yellow are the 95 % HPD. A time scale show the geological ages. Gray shaded areas represent the GPG (~1.68-1.016 Mya) and the CPG (~0.7-0.6 Mya). Below the tree glacial-interglacial cycles were reconstructed using the data available from [82] and the historical ocean temperature available from [84]

Discussion

Geographical distribution of Clades

We recovered four groups at the north of the Maipo river (southern clade): (1) The Loa lineage, (2) T. e. coquimbensis, (3) T. e. elegans and (4) the new lineage Aconcagua. The Loa lineage (Clade A) is located north of the Loa river mouth (21°25’S), and so far is the only population known for this lineage, but specimens have been recorded in Paposo-Taltal [30] which could belong to this lineage as well. Whereas T. e. coquimbensis (Clade B) has the widest distributional range. Our results extend its distributional range between the coast of Atacama southward to the Coquimbo region from Quebrada El León (26°57’S) to Valle El Mauro (31°58’S), with the Quilimarí river and Santa Inés-Imán mountain range as southern boundary (Figs. 2 and 3, Fig. S1 and Fig. S2). These two lineages inhabits in the hyperarid Atacama in localities as the northern Loa river mouth, Quebrada El León and Llanos del Challe where the rain barely reaches the 10 to 20 mm/yr [75, 76] and are the most extremophilous lineage which inhabit in small patches, strongly associated to rivers, ravines, or riverbeds taking advantage of the available water and the small bushes present there (see Fig. 1a, b).

On the other hand, Our results showed that in the next southern clades the ranges of distribution are drastically reduced. The new lineage Aconcagua (Clade C) is restricted from Fundo El Roble (32°16’S) to Fundo Chuico Blanco (32°48’S), between the Quilimarí river-Santa Inés-Imán mountain range to the Aconcagua river and the Chacabuco mountain range. Same than T. e. elegans (Clade D) which is restricted between the Aconcagua river and Chacabuco range (32°57’S) and the Maipo river and Altos de Cantillana range (33°34’S). At the South of the Maipo river (southern clade), two lineages previously not described were recovered. The South 1 lineage (Clade E) between the Maipo river and Altos de Cantillana range (33°3’S) and the Mataquito river (35°3’S) and its tributary Teno river (34°55’S). Finally, the South 2 lineage (Clade F), restricted at south of the Mataquito-Teno rivers at least until Tregualemu (35°56’S). Interestingly, the specimens from the South 2 lineage could be assigned to T. e. soricinus according to [29], who describe its distribution until Radal Siete Tazas and the junction of the Claro and Maule rivers, which is totally consistent with the distribution of the South 2 lineage (Figs. 2 and 3, Fig. S1 and Fig. S2). Historical records [21, 22, 32] propose that T. e. soricinus could inhabit as far as the south of the Bío-Bío river until Angol and adjacent areas of the Nahuelbuta Coastal range 37°50’S (see Fig. 1f).

In the northern part of the Mediterranean region (31°–35°S), thorny scrublands alternate with sclerophyllous woodlands and relict populations of the Chilean palm, in the Valparaíso region, while in the southernmost zone (36°–38°S), sclerophyllous forests mixed with broadleaf evergreen and deciduous forests form a diverse ecotone in the Mediterranean-temperate climatic transition [26]. In summary, from the landscape point of view, we can divide the distribution of T. elegans lineages in three climatic zones adapted from [26,27,28]: (1) Desert characterized by scarce or null vegetation where inhabit mainly the clades A in Tarapacá region and northern part of clade B in Atacama region, (2) Dry Mediterranean with vegetation type such as dry xerophytic thorn scrublands and evergreen sclerophyllous communities where inhabit mainly the clades C, D and E, but also the southern populations of clade B in Coquimbo region and 3) Wet Mediterranean, dominated by deciduous forest were inhabiting the clade F in Maule region.

Phylogeographic determinants

About 3.5 Mya a dramatic climatic change occurred because was the end of the warm Pliocene as result of global cooling. During this major climate cooling, important glaciations in all Andes occurred [77, 78] and consequently the rivers formation post-Pliocene glaciation [79, 80]. Following further cooling towards the Quaternary and throughout the Quaternary established the South American Arid Diagonal [81] which change the landscape of the continent as is known today. Likely, all these climatic events triggered the separation of the northern and southern Maipo clades by a vicariant event, because the Maipo river flow increase posterior to a glaciation and cooling periods due to melting of glaciers. These geographic barriers could be reinforced for the transversal range that joined the Andes with the Coastal Cantillana range, the Maipo river has been a hard barrier from this period. The reciprocally monophyletic clades in opposite banks and the oldest divergence time (3.14 Mya) are strongly supporting the primary vicariant event as hypothesis of diversification on the Pleistocene (Figs. 2 and 4, Fig. S1, Fig. S2 and Fig. S3). After a huge gap of diversification, around 1.8 Mya stands for the most fundamental change in global climate and climate variability in the Quaternary, referred to as the Mid Pleistocene Transition [82,83,84]. This period coincides with the significant emergence of diversification events in T. elegans at ~1.8-0.8 Mya. During this period, the previous dominant periodicity of climate cycles changed from a 41 kyr to a 100 kyr cycle, causing high amplitude climate oscillations (e.g.[82, 85]). Pronounced cold and warm cycles during the remaining Quaternary most likely caused the repeated cyclical opening and closing of migration corridors, in this case by rivers or hyperaridity, promoting further differentiation in the lineages of T. elegans. These divergence time match with the deep effect that the Great Patagonian Glaciation (GPG~ 1.68-1.016 Mya) generated on a wide range of taxa in southern South America [86]. Therefore, glacial-interglacial cycles of the Quaternary could affect the river water levels fluctuations of the Mediterranean Chilean zone, becoming potential barriers to gene flow of the lineages of T. elegans during melting of glaciers and more permeable barriers during glaciations[35, 36, 40, 55, 87, 88]. During glacial- interglacial cycles, great ice masses descending to around 1100-1300 masl between 33-34°S in the Aconcagua and Maipo valleys [89,90,91]. Around ∼1.84 Mya, the T. e. elegans lineage was separated from the northern lineages because the Aconcagua river and Chacabuco range acted as effective geographic barriers and biogeographical determinants to keep its ancestral range. While the other lineages of T. elegans differed in central and northern Chile in a relatively short period of time. Around ∼1.41 Mya, the lineage T. e. coquimbensis occupied northern areas probably because a hyperarid period between 2 and 1 Mya affected the Atacama Desert region [92], allowing that a small population crossed the Quilimarí river though the coast because the lower river water levels in middle courses and river mouths of the Mediterranean Chilean zone were more permeable barriers [93].

Then, the peripherical differentiation in the Coquimbean Andes avoiding the hyperaridity of the desert, facilitated that this lineage reach its northernmost distribution along the coast in Quebrada el León in posterior favorable conditions. Today the Quilimarí river and the Santa Inés-Imán mountain range seem to be the southern distribution boundary of T. e. coquimbensis but it is not clear if they are strong barriers to southern dispersion allowing secondary divergence. Finally, one of the most complex scenarios occurred about 1.28 Mya, where Central Chile was colonized. Some populations from the Atacama Desert at the north of the Loa river mouth, differentiated into the Loa lineage. Likely, the Loa lineage had a distribution extended further south in the Atacama Desert. Then, the populations were connected by old corridors and perhaps the southernmost dispersal to Central Chile. Huge corridors in the Atacama Desert have been proposed for other organism to explain this strange vicariant distribution with a huge gap [94, 95]. However, these corridor are not permanent and probably have been severely closed during hyperarid periods (as it is now) in the Atacama and today the biota can survive in small pocket with suitable conditions [94], while other middle population were extinguished. The Aconcagua lineage finally differentiated between the Quilimarí river and Santa Inés-Imán mountain range and the Aconcagua river and the Chacabuco range as barrier to southern dispersion allowing secondary divergence. Probably, these events are related, occurring almost at the same period.

On the other hand, the southern history is more recent and simpler, because the phylogeny recovered reciprocally monophyletic groups on opposite banks suggesting that around 0.81 Mya during the Pleistocene, a vicariant event promoted by the Mataquito river and its tributary Teno river caused the divergence between the South 1 and South 2 lineages (Figs. 2 and 4, Fig. S1, Fig. S2 and Fig. S3). Chile was strongly affected by the Pleistocene glaciation, particularly in the south part of the country with an event as the Great Patagonian Glaciation (GPG~1.68-1.016 Mya) and the Coldest Pleistocene Glaciation (CPG~0.7-0.6 Mya). These glaciations could be related with the water levels increase of the Mataquito and Teno rivers [96] and acting such as primary vicariant barrier between both lineage. Mataquito river has been reported as a phylogeographic break at 35° S for terrestrial animals [38, 42].

In general, we have two spatially contrasting scenarios for explain the structure found. Central Chile understanding as the region from the Quilimarí river (32°6’S) to the Bío-Bío river (36°48’S) had a greater influence of glaciations, generating a greater water levels in the past for rivers [35, 36, 40, 55, 87, 88] and these rivers have strongly influenced on the genetic structure of T. elegans, but with different intensity. Thus, Maipo and Mataquito-Teno rivers acted such as primary and strong barriers, while rivers north to Maipo such as Quilimarí and Aconcagua are secondaries barrier which in some parts of their history were permeable for the fauna. On this way, the past climate changes to the hyperaridity would be the key to explain the current distributional range of the T. elegans lineages in the Atacama Desert and northern Mediterranean zone of Chile (i.e. north of the Quilimarí river). The Loa lineage which is enclosed in a small pocket in the north Loa river mouth is currently distributed far of its sister clade (Aconcagua lineage) in Central Chile and is the good example of the dramatic effect of the aridity for the biota, probably extinguishing populations between the Loa and Aconcagua lineages (Fig. 2, Fig. S1 and Fig. S2). However, T. e. coquimbesis would have recently replaced whose anteriorly extinguished populations from Coquimbean Andes. The huge distribution of T. e. coquimbesis in the Atacama desert (Figs. 2 and 3), without genetic structure seem indicate that the effects of rivers are apparently null there, even when some rivers have been reported as important barrier for animals, particularly the Huasco and Copiapó rivers [51, 55]. However, the survivor Atacama populations today are living in small patches with vegetation which are small island in the hyperarid Atacama.

Probably with the creation of the Arid Diagonal and the periodic climate oscillation which changed the landscape several times in the last 1.5 Mya, the different ancestors of each lineage were adapted to xeric environment (clades A and northern part of clade B), mesic (clades C, D,E and southern part of clade B) or humid environment (Clade F) favoring even more the differentiation not drove by barriers. A similar scenario was reported in Abrothrix longipilis with lineages restricted to different bioclimate [26, 47]. Therefore, although today some rivers as the Quilimarí or even the Aconcagua seems to be barriers, that they could be only delimitating the current distribution of the lineages but not promoted divergence by vicariance as suggest the topology and phylogenetic relationships of our results because T. e. coquimbensis and Aconcagua lineage from opposite banks of Quilimarí and Aconcagua lineage and T. e. elegans from opposite banks of Aconcagua river are not sister clades (Fig. 2, Fig. S1 and Fig. S2). Probably, the current population of each clade is enclosed in particular patches avoiding the aridity. Therefore, the scenario of permeable barriers is also possible, but have not been detected due to the elapsed time or because these lineages have not crossed the rivers since they are in its habitat distribution avoiding the aridity. In the future its necessary to evaluate explicitly the synergistic effect of topography, glaciations and rivers in the phylogeographic structure of more Chilean species from Mediterranean zone of Chile. The Andean and Coastal range, multiple ice sheet glacial advances and retreats and increases in level water volumes due to melting ice probably influenced terrestrial communities in this region. As such, widely distributed endemic species as T. elegans is a good model to assess the joint effects of topography, glacial cycles, and changes in the river water volumes on its genetic differentiation and phylogeographic structure of T. elegans lineages [36, 52, 97].

Taxonomic implications

Apparently, the intraspecific diversity of T. elegans has been historically underrated, and our evidence shows six clearly differentiated genetic lineages. The Loa lineage was previously reported [19, 20], meanwhile a drastic reduction of distributional range was reported in T. e. elegans and an increase to coastal areas from Atacama and Coquimbo region in T. e. coquimbensis [20]. However, the Aconcagua, South 1 and South 2 are newly lineages reported in this study (Figs. 2 and 3, Fig. S1 and Fig. S2).

In general, the genetic distance for species of the genus Thylamys are between 2.5 and 20 % [18, 19, 98]. In our case, the values observed in the genetic differences for the six lineages of T. elegans was between 2 and 11.9 % For example, we are reporting a genetic distance of 9.7 % between the northern and southern clades (Table 1) which is even higher than 5.4 % between T. elegans “north” and T. elegans “south” detected by [19] or higher than linages found within of T. pallidior (5 %), T. venustus (5.4-3.9 %), and T. sponsorious (2.5 %) [18]. However, a study with 82 marsupial species reported that 95 % of cyt b interspecies genetic distances in marsupials have values between 9 % and 16 % with a mean of 12.3 % [99].

Therefore, our results are suggesting a cryptic taxonomic scenario which should be approached by integrative evidence considering species delimitation species analyses, other genes and/or morphological characters. The evidence presented is only the starting point to highlight the deep genetic structure within T. elegans and for the moment, our study is a simple call for futures accurate taxonomic studies on this endemic mouse opossum.

Conclusions

Biogeographic history of T. elegans was shaped mainly for changes of the rivers water levels during the Pliocene and Pleistocene glaciations modulating the genetic differentiation of its lineages. Our results sustain a comprehensive overview suggesting how and when the genetic diversity in T. elegans was structured in the past. We found different intensity of the phylogeographic determinants in the northern and southern clades. For the first time we reported the effects of the Maipo, Mataquito and Teno rivers on the genetic structure and lineage divergence of an endemic small mammal of Chile as primary riverine divergence. While the Quilimarí and Aconcagua rivers were barriers to dispersal, acting as actual geographic boundary and geographic barriers preventing secondary contact in periods of demography and geographic expansion such a secondary riverine divergence. Therefore, rivers must be considered as dynamic barriers that fluctuate over time, have been effective barriers to gene flow. However, the Maipo river is a strong barrier which produced the deep divergence between the northern and southern clades.

Other important barriers could be the final uplift of the Andes, the mountain ranges and the transverse valleys formation. We also found a synergistic role of the climatic changes, hyperaridity, glaciations and mountain ranges formations as biogeographic factors, in the origin and preservation of the genetic divergence of T. elegans suggesting effects of several events and a complex regional history in Chilean Mediterranean zone. All these factors structured deeply the genetic diversity into T. elegans which is composed by six lineages distributed latitudinally and whose genetic distance is suggesting a cryptic complex. This study is only an open window for deeper and detailed studies into this endemic mouse opossum.