The riverine barrier model suggests that rivers play a significant role in separating widespread organisms into isolated populations. In this study, we used a comparative approach to investigate the phylogeography of 6 didelphid marsupial species in central Brazil. Specifically, we evaluate the role of the mid-Araguaia River in differentiating populations and estimate divergence time among lineages to assess the timing of differentiation of these species, using mitochondrial DNA sequence data. The 6 didelphid marsupials revealed different intraspecific genetic patterns and structure. The 3 larger and more generalist species, Didelphis albiventris, Didelphis marsupialis, and Philander opossum, showed connectivity across the Araguaia River. In contrast the genetic structure of the 3 smaller and specialist species, Gracilinanus agilis, Marmosa (Marmosa) murina, and Marmosa (Micoureus) demerarae was shaped by the mid-Araguaia. Moreover, the split of eastern and western bank populations of the 2 latter species is consistent with the age of Araguaia River sediments formation. We hypothesize that the role of the Araguaia as a riverine barrier is linked to the level of ecological specialization among the 6 didelphid species and differences in their ability to cross rivers or disperse through the associated habitat types.

The riverine barrier model was first proposed by Wallace (1852), suggesting that rivers may have a significant role in separating widespread organisms into isolated populations. Several vertebrate groups were extensively sampled in a pioneer study in the Juruá River in Amazonia aiming to examine the riverine effects on patterns of both population and species differentiation and ecological distribution (Patton et al. 2000). This and several other successive studies that investigated the potential role of Neotropical rivers as putative barriers reached different conclusions (e.g., Ayres and Clutton-Brock 1992; Gascon et al. 2000; Matocq et al. 2000; Bates et al. 2004; Faria et al. 2013a; Maldonado-Coelho et al. 2013). The lack of congruence among species and rivers has been related with the species vagility and also with the drainage system geography and formation history (Moritz et al. 2000; Antonelli et al. 2010; Leite and Rogers 2013).

The Araguaia River lies in the transitional area between Cerrado and Amazonia. With a length of 2110 km and reaching 3–6 km of width in its middle portion, it represents the largest river of the wet-dry tropics of Brazil (Latrubesse and Stevaux 2002). This river is classified as an anabranching river, consisting of channels separated by vegetated semi-permanent alluvial islands, with low sinuosities (Latrubesse 2008). Like other rivers in the Brazilian Shield it represents a relatively stationary geographic feature, which is expected to limit gene flow between forest-dwelling populations from opposite banks, contrary to meandering rivers in the western region of the Amazon Basin (e.g., Aleixo 2006; Bates et al. 2004; Antonelli et al. 2010; Maldonado-Coelho et al. 2013). Therefore, the Araguaia River provides an excellent natural setting to test biogeographic diversification hypotheses such as the riverine barrier model. Previous studies reached different conclusions, showing that Araguaia might have been a barrier to 2 Rhipidomys rodent sister species (Rocha et al. 2011), but not playing an important role in the intraspecific differentiation of the rodents Hylaeamys megacephalus and Oecomys aff. roberti (Rocha et al. 2014). Thus, it is essential to understand to what extent these findings are consistent across different animal taxa, that is, in a general biogeographic framework, and why different species are responding differently to the same geographic feature.

The family Didelphidae is among the oldest of extant mammal families (Jansa et al. 2014). Species of this diverse group are common in forest and open habitat communities in the Neotropics (Gardner 2008). Currently, 11 species are recognized in the ecotonal area between Cerrado and Amazonia (Lacher and Alho 2001; Bezerra et al. 2009; Rocha et al. 2011). In this study, we considered the 6 most common species of didelphid marsupials observed in the mid-Araguaia River basin (Rocha et al. 2011): Didelphis albiventris, Didelphis marsupialis, Philander opossum, Gracilinanus agilis, Marmosa (Marmosa) murina and Marmosa (Micoureus) demerarae (for recent taxonomic revision see Voss et al. 2014). Although these 6 species are widely distributed in lowland Neotropical forests (Gardner 2008), few studies have focused on their phylogeography in central Brazil (exceptions include Costa 2003; Faria et al. 2013a, 2013b), and none in the ecotonal area between Cerrado and Amazonia, which remains poorly studied overall (but see Lacher and Alho 2001; Bezerra et al. 2009; Rocha et al. 2011).

Here, mitochondrial DNA (mtDNA) sequence data was used to investigate the role of the Araguaia River in differentiating populations of 6 didelphid marsupials in central Brazil. Inferring the evolutionary and demographic history of species solely on mtDNA can be misleading, due to its strictly maternal inheritance (Ballard and Whitlock 2004), but exploratory analyses with this marker are important to provide basic knowledge of putative barriers and geomorphological events (e.g., Patton et al. 2000; Costa 2003). Furthermore, molecular data combined with geographic distribution information have been useful to understand the evolutionary history of species and to explicitly test hypotheses of biogeographical events (e.g., Patton et al. 2000; Costa 2003; Antonelli et al. 2010; Nicolas et al. 2011). Additionally, comparative phylogeography, using multiple codistributed taxa, allows a better understanding of the relationships between landscape formation and biotic diversification (Arbogast and Kenagy 2001).

The distribution of haplotypes across the river holds information about the population processes that have shaped the history of each species (e.g., Patton et al. 2000; Nicolas et al. 2011). Patton et al. (2000) formulated 3 different hypotheses for riverine divergence, which included primary diversification, secondary contact, and dispersal. The first hypothesis postulates a complete barrier formed by the river to an existing taxon range, resulting in a reciprocally monophyly of sister clades from opposite banks. In the secondary contact hypothesis, the river is only a common contact zone, and despite of the observed monophyly, haplotypes from opposite banks are not sister clades. The third case comprises dispersal events from an established population to the opposite riverbank, resulting in paraphyletic relationships. Moreover, low differentiation will be revealed by a certain degree of haplotype sharing across the river (see Patton et al. 2000 for further explanations). To frame the phylogenetic predictions of the river barrier hypothesis we included in the analysis samples of didelphid marsupial species from other geographic areas, including the Atlantic Forest, Pantanal, and Amazonia. We also estimated divergence time among lineages to assess the timing of differentiation and to assess the potential impact of external events (e.g., geomorphological) on population differentiation (Leite and Rogers 2013). If the mid-Araguaia River represents a barrier to populations of didelphid species, it is expected that the diversification between populations from opposite river banks coincides with the formation of the river during the Pleistocene, roughly 240000 years ago (Latrubesse and Stevaux 2002; Valente and Latrubesse 2012). This means that for older divergence events, the river is unlikely to be the primary agent of divergence, although it could nonetheless be acting as barrier to gene flow (e.g., Maldonado-Coelho et al. 2013). Finally, we also explored the species demographic history to detect signatures of population expansion and/or bottlenecks. If the mid-Araguaia River was a primary barrier in population divergence, one should not find evidence of population expansion by the river (Moritz et al. 2000; Leite and Rogers 2013). On the other hand, if there is signature of historical population fluctuations, the Araguaia River could be just a secondary barrier as populations could have been isolated by other geographic processes (i.e., refuges; Moritz et al. 2000; Leite and Rogers 2013).

Methods

Sample Collection

This work focused on 6 didelphid marsupials: D. albiventris (n = 15), D. marsupialis (n = 50), G. agilis (n = 108), M. murina (n = 80), M. demerarae (n = 35), and P. opossum (n = 54) from central Brazil. Samples were collected in several localities in the both margins of the mid-Araguaia River, 2 islands within this river and also in its tributaries, Javaés and Coco Rivers (Figure 1). Sampling was carried out between June 2007 and November 2008, using a standardized trapping protocol to sample small nonvolant mammals in upland and floodplain gallery forests (for detailed sampling at mid-Araguaia see Rocha et al. 2011). Tissue samples of the mid-Araguaia River have been deposited at Universidade Federal do Espírito Santo (UFES), Vitória, Brazil. Additional samples from central Brazil were obtained via museum tissue collection or kindly shared by colleagues (for details on biological material and sequences used, including sampling localities for each species see Supplementary Table S1 and Gazetteer online). The central Brazilian region treated in this work encompasses mainly samples from the ecotonal region between Cerrado and Amazonia, but also samples from Cerrado (Figure 1).

Figure 1.

(a) Sampling areas at the mid-Araguaia River and its tributaries (left, source: Google Earth). Sampling points in white corresponds to the eastern river bank, in black to the western river bank, and in gray to the islands of the Araguaia River. (b) Map including all sampling localities of 6 didelphid marsupials in South America (right). See Gazetteer for detailed information on localities.

Figure 1.

(a) Sampling areas at the mid-Araguaia River and its tributaries (left, source: Google Earth). Sampling points in white corresponds to the eastern river bank, in black to the western river bank, and in gray to the islands of the Araguaia River. (b) Map including all sampling localities of 6 didelphid marsupials in South America (right). See Gazetteer for detailed information on localities.

Additionally, we included sequences of didelphid marsupial species from other geographic areas, including the Atlantic Forest, Pantanal, and Amazonia (Figure 1), available in GenBank: D. albiventris (n = 5), D. marsupialis (n = 8), G. agilis (n = 26), M. murina (n = 23), M. demerarae (n = 20), and P. opossum (n = 6), in order to investigate the geographic structure and phylogeographic affinities of central Brazil.

DNA Extraction and Amplification

DNA was extracted from liver or ear tissue preserved in ethanol using salt-extraction method (Bruford et al. 1992). Cytochrome b (mt-Cytb) fragments with 801 base pairs (bp) were amplified by polymerase chain reaction (PCR) using the primers MVZ05 and MVZ16 (Smith and Patton 1993). PCR reactions were performed in a 25 μL total volume, including 2.5 μL of 10× PCR buffer, 1.0 μL of MgCl2 (50mM), 0.5 μL of dNTP mixture (10mM for each nucleotide), 0.3 μL of Platinum® Taq DNA polymerase (Invitrogen, Waltham, MA)), 0.3 μL of each primer (10 μM), and 1 μL of DNA template (100ng/μL). Amplifications were performed with the following PCR profile: initial denaturation at 94 °C (5min), followed by 39 cycles with denaturation at 94 °C (30 s), annealing at 48 °C (45 s), polymerization at 72 °C (45 s), and a final extension at 72 °C (5min). Mitochondrial fragments were purified using ExoSap-IT® (USB Corporation) and sequenced using an automatic sequencer ABI 3130-XL (Perkin Elmer, Applied Biosystems, Foster City, CA), with the above-listed primers. Sequence alignment was performed using CLUSTALW algorithm implemented in MEGA 6 (Tamura et al. 2013). Electropherograms and sequences were manually checked and edited in MEGA 6.

Data Archiving

In accordance with the data archiving policy (Baker 2013), all sequences generated in this study have been deposited in GenBank (see Supplementary Table S1 online).

Phylogeographic Relationships

In order to investigate the geographic distribution and phylogeographic affinities of 6 didelphid marsupials, and specifically to test if the mid-Araguaia River has a significant role in structuring populations of those didelphid species in central Brazil, relationships amongst haplotypes were estimated through Bayesian inference (BI) performed in the MrBayes v.3.1.2 (Ronquist and Huelsenbeck 2003). For this analysis, we used only one individual per haplotype representing populations from central Brazil, but also haplotypes from the Atlantic Forest, Pantanal, and other regions in the Amazonia were included in posterior analyses. Haplotypes from closely related species were used as outgroups. The best model of nucleotide substitution was selected in MrModeltest (Nylander 2004). Haplotype trees were sampled every 100 of 107 generations until Markov chain became stationary, that is, when standard deviation of split frequencies was below 0.01. A 50% majority rule consensus tree was obtained after burn-in of 25% of the sample points to generate Bayesian posterior probabilities (BPP).

Pairwise sequence distances (%) within and between the main geographic clades were estimated using the Kimura 2-parameter (K2p) model (Kimura 1980) implemented in MEGA 6.

Median-joining (MJ) networks were constructed for central Brazilian samples only (see Supplementary Table S1 online), using NETWORK v.4.6 (Bandelt et al. 1999). Only variable nucleotide sites were included in MJ analyses. For the purposes of this analysis, we categorized 2 areas as comprising all populations from west and east of the Araguaia River. Haplotypes from each bank of the Araguaia River were identified with different colours (white—eastern bank, EAR; black—western bank, WAR) to illustrate the riverine barrier or put in evidence the haplotype sharing across the river. For D. marsupialis, samples from islands of the Araguaia River were also available and we categorized a third area for this species, identifying samples from those islands with gray colour.

Divergence Time and Gene Flow

To estimate the time to most recent common ancestor (tMRCA) we ran a single data set with 87 mt-Cytb sequences including haplotypes from all 6 species from this study (see Supplementary Table S1 online) and, in order to use multiple calibration points, we included other Didelphidae sequences obtained from GenBank (Caluromys lanatus: U34664; Chironectes minimus: AJ628363; Lutreolina crassicaudata: AJ628364; Monodelphis domestica: HQ651773; Thylamys karimii: EF051700). Molecular clock was tested by performing a maximum likelihood (ML) analyses in MEGA 6. The ML values of the obtained topologies were compared with and without molecular clock constraints also using MEGA 6. The presence of strict molecular clock for this dataset was rejected. The tMRCA analyses were performed in BEAST 2.1.3 (Bouckaert et al. 2014) using the Bayesian relaxed clock model and allowing the branch length to vary according to an uncorrelated lognormal distribution (Drummond et al. 2006). We used the Yule process speciation model as tree prior and the GTR+I+G evolution model. We followed the phylogeny of Voss and Jansa (2009) to set monophyletic constrains. Four calibration points were used: the diversification of the Didelphidae and the diversification of Didelphimorphia, which were set as minimum and maximum ages (Meredith et al. 2011), and 2 minimum divergence dates within the marsupial radiation, which were set assuming a log-normal distribution with standard deviation of 1.0 (Jansa et al. 2014). The minimum age of the diversification of Didelphidae and Didelphimorphia was the oldest unequivocal fossil belonging to clade at 11.608 million years ago before present (Myr BP) for both (Meredith et al. 2011). The maximum age of Didelphidae was based on stratigraphic bounding, which encompasses 2 stratigraphic units without fossil records of this clade at 28.5 Myr BP (Meredith et al. 2011). The maximum age of Didelphimorphia was based on phylogenetic bounding, which encompasses the age of the stem of didelphimorph at 65.8 Myr BP (Peradectes) (Horovitz et al. 2009; Meredith et al. 2011). The 2 minimum divergence dates within the marsupial clade were the divergence between Didelphis and Philander at 3.3 Myr BP; and the split between Monodelphis and Marmosa at 12.1 Myr BP (Jansa et al. 2014). We performed 2 independent runs of 100000000 generations each and sampled every 10000 generations. Analyses were performed in the CIPRES Science Gateway (https://www.phylo.org). Convergence of the MCMC chains was verified in Tracer v.1.6 by checking the effective sample size values (ESS). The resulting trees were combined using LogCombiner and TreeAnotator after a burn-in of 10%. Consensus trees were visualized in FigTree 1.4.

The IM program (Hey and Nielsen 2004) was used to estimate the effective number of migrants per generation (m), time of divergence (t) between populations from either river bank of the Araguaia River, and the effective population size in both populations as well as the ancestral population (qEAR, qWAR, and qA, respectively). For these analyses we only used samples from the immediate vicinity of the Araguaia River (see Supplementary Table S1 online). For G. agilis we did not perform this analysis since we had only a single sample from the western river bank. For D. marsupialis we categorized 3 populations (as this species inhabits both the river banks and mid-river islands) and we performed all 3 pairwise combinations: eastern and western river banks, eastern river bank and islands and western river bank and islands. We used an inheritance scale of 0.25 for the mtDNA and the HKY model of evolution. Substitution rates per nucleotide were obtained using BEAST (as mentioned above) and transformed into complete-locus substitutions, resulting in mean of 0.000014 and 95% highest posterior density (95% HPD) of 0.000003–0.00003. We assumed a generation time of 1 year for all species. A single unheated chain was used to explore the parameter space for each species. We did at least 5 preliminary trial runs for each species in order to assess the range of prior distributions prior for theta (varying from 10 to 30), migration (varying from 2 to 20) and time estimation (varying from 6 to 20) (Hey and Nielsen 2004). We assumed symmetrical gene flow (i.e., m1 = m2) for all species to avoid over-parameterizing the model. We used a burn-in period of 500000 steps and recorded results every hour. We allowed the IM program to run until the lowest ESS were at least 500 (e.g., Peters et al. 2005). We visually inspected the parameter trend lines over each MCMC chain and recorded ESS of the recorded values for each parameter. We only reported results from the longest run for all species (D. albiventris: 64217790 steps, minimum ESS = 5570; D. marsupialis: 10663380 steps, minimum ESS = 4548; P. opossum: 8258956 steps, minimum ESS = 5027; M. murina: 32290531 steps, minimum ESS = 1388; and M. demerarae: 51902032 steps, minimum ESS = 14564).

Demographic History

Number of haplotypes (h), number of polymorphic sites (S), haplotype (Hd) and nucleotide (π) diversity values for the central Brazilian samples were estimated with DnaSP v.5 (Librado and Rozas 2009).

Deviation from neutrality in central Brazil was tested through the following statistics based on frequency spectrum of mutations, Tajima’s D (Tajima 1989) and Ramos-Onsins and Rozas R2 (Ramos-Onsins and Rozas 2002), and based on haplotype frequencies, Fu’s Fs (Fu 1997), using DnaSP v.5. Coalescence simulations with 1000 replicates were applied to determine the P-value of each statistics. Significant P-values (<0.05) were taken as evidence a scenario of demographic expansion.

Demographic history of the 6 didelphid marsupials was further investigated using Bayesian Skylines Plots (BSP) implemented in BEAST 2.1.3. We used substitution rates per nucleotide per million year obtained in the above-mentioned estimates of the divergence time (Mean = 0.017, 95% HPD = 0.004–0.035), and prior best models of nucleotide substitution selected in MrModeltest. Two independent runs of 100000000 generations each and sampled every 10000 generations, were performed for each analysis. Convergence of the MCMC chains was verified in Tracer v.1.6 by checking the ESS values. Skyline plots were constructed using Tracer v.1.6.

Results

Didelphis albiventris

The BI analyses based on 14 haplotypes showed 3 well-supported clades for D. albiventris: the Southwest clade (SWC), the Cerrado clade (CE), and the Amazonia-Cerrado ecotone clade (ECO) (Figure 2a). The SWC clade included haplotypes from Southern Brazilian Atlantic Forest and from Bolivian Chaco domains (Figure 2). The ECO clade is represented mainly by haplotypes from the Amazonia-Cerrado ecotone, especially from mid-Araguaia River (Figure 2). The CE clade is represented mainly by haplotypes from Cerrado (Figure 2). However, haplotype HDa1 from the Amazonia-Cerrado ecotone (locality 34) is closely related to the CE clade, and haplotype HDa8 from Cerrado (locality 77) is closely related to the ECO clade (Figure 2).

Figure 2.

(a) Bayesian phylogenetic reconstruction of mt-Cytb haplotypes of Didelphis albiventris (GTR+I model used). Numbers in circles represent nodes for which time to most recent common ancestor are shown in Table 1; percentage values correspond to pairwise genetic distances between clades; and asterisks indicate BPP ≥0.95; numbers in the tip of branches are sampling localities and acronyms are haplotype designations; clades are designated by ECO (Amazonia-Cerrado ecotone clade), CE (Cerrado clade), and SWC (southwest clade). (b) MJ network of D. albiventris of the ECO and CE clades, central Brazil. Length of connecting branches corresponds to the nucleotide substitutions; size of the circles is proportional to the number of individuals sharing each haplotype; white corresponds to the eastern riverbank of the Araguaia River, and black to the western riverbank. (c) Map of sampling localities of D. albiventris, with the representation of the clades.

Figure 2.

(a) Bayesian phylogenetic reconstruction of mt-Cytb haplotypes of Didelphis albiventris (GTR+I model used). Numbers in circles represent nodes for which time to most recent common ancestor are shown in Table 1; percentage values correspond to pairwise genetic distances between clades; and asterisks indicate BPP ≥0.95; numbers in the tip of branches are sampling localities and acronyms are haplotype designations; clades are designated by ECO (Amazonia-Cerrado ecotone clade), CE (Cerrado clade), and SWC (southwest clade). (b) MJ network of D. albiventris of the ECO and CE clades, central Brazil. Length of connecting branches corresponds to the nucleotide substitutions; size of the circles is proportional to the number of individuals sharing each haplotype; white corresponds to the eastern riverbank of the Araguaia River, and black to the western riverbank. (c) Map of sampling localities of D. albiventris, with the representation of the clades.

The average genetic distance between SWC clade and CE plus ECO clades is 6.6±0.9%, and their tMRCA dates to the Pleistocene at 0.76–2.52 Myr BP (Table 1, Supplementary Figure S1 online). The average K2p genetic distance between the 2 central Brazilian clades (CE and ECO) is of 1.1±0.3%, and their tMRCA dates to the Pleistocene at 0.24–0.98 Myr BP (Table 1, Supplementary Figure S1 online).

Table 1.

Mean time to most recent common ancestor of the main clades of didelphid species

Node Mean 95% CI Geological time 
1 (Didelphis albiventris ECO / D. albiventris CE) 0.56 0.24–0.98 Pleistocene 
2 (D. albiventris ECO + CE / D. albiventris SWC) 1.56 0.76–2.52 Pleistocene 
1 (Didelphis marsupialis0.98 0.44–1.65 Pleistocene 
1 (Philander opossum1.11 0.49–1.87 Pleistocene 
1 (Gracilinanus agilis WCE + CBP + CEC + SC) 1.26 0.71–1.94 Pleistocene 
2 (G. agilis WCE + CBP + CEC + SC / G. agilis NCE) 1.50 0.83–2.32 Pleistocene 
3 (G. agilis2.15 1.15–3.32 Pliocene/Pleistocene 
1 (Marmosa murina EAR / M. murina WAR) 0.81 0.39–1.32 Pleistocene 
2 (M. murina1.99 1.00–3.19 Pliocene/Pleistocene 
1 (Marmosa demerarae ECO + EToR + AF) 1.35 0.69–2.10 Pleistocene 
2 (M. demerarae ECO + EToR + AF / M. demerarae GS + NAM) 1.76 0.96–2.68 Pliocene/Pleistocene 
3 (M. demerarae2.36 1.34–3.61 Pliocene/Pleistocene 
Node Mean 95% CI Geological time 
1 (Didelphis albiventris ECO / D. albiventris CE) 0.56 0.24–0.98 Pleistocene 
2 (D. albiventris ECO + CE / D. albiventris SWC) 1.56 0.76–2.52 Pleistocene 
1 (Didelphis marsupialis0.98 0.44–1.65 Pleistocene 
1 (Philander opossum1.11 0.49–1.87 Pleistocene 
1 (Gracilinanus agilis WCE + CBP + CEC + SC) 1.26 0.71–1.94 Pleistocene 
2 (G. agilis WCE + CBP + CEC + SC / G. agilis NCE) 1.50 0.83–2.32 Pleistocene 
3 (G. agilis2.15 1.15–3.32 Pliocene/Pleistocene 
1 (Marmosa murina EAR / M. murina WAR) 0.81 0.39–1.32 Pleistocene 
2 (M. murina1.99 1.00–3.19 Pliocene/Pleistocene 
1 (Marmosa demerarae ECO + EToR + AF) 1.35 0.69–2.10 Pleistocene 
2 (M. demerarae ECO + EToR + AF / M. demerarae GS + NAM) 1.76 0.96–2.68 Pliocene/Pleistocene 
3 (M. demerarae2.36 1.34–3.61 Pliocene/Pleistocene 

Mean time and 95% confidence intervals (95% CI) are given in Myr BP. Geological time correspondent to the tMRCA adapted from Cohen et al. (2013).

The MJ network was constructed with 15 samples from central Brazil (Table 2). CE and ECO clade are separated by 4 mutations (Figure 2b). Haplotype HDa4 is the ancestor of the CE clade, and it is shared by specimens from both riverbanks of the Araguaia (Figure 2b). Hypothetical ancestor of the ECO clade was not sampled (Figure 2b). The migration rate m estimated in IM (using ECO clade) had a long tail toward higher values (see Supplementary Figure S2 online). The best estimate of divergence time t between the eastern and western river banks was 0.645 Myr BP (95% HPD = 0.031–1.389 Myr BP, Table 3). Eastern (qEAR), western (qWAR), and ancestral (qA) populations had similar inferred sizes (Table 3, Supplementary Figure S3 online). Neutrality tests were not significant (Table 2), that is, they did not reject the constant population model. BSP for the ECO clade also support a constant population (Figure 3).

Table 2.

Number of individual sequences (n), haplotypes (h), number of polymorphic sites (S), haplotype (Hd), and nucleotide (π) diversity, and deviation from neutrality tests (Tajima’s D, Fu’s Fs, and Ramos-Onsins and Rozas’s R2) of 6 didelphid marsupial cyt-b sequences for central Brazil

Species n h S Hd π D FR2 
Didelphis albiventris 15 18 0.876 0.0062 −0.40 NS −1.112NS 0.1246NS 
Didelphis marsupialis 50 0.642 0.0014 −0.378NS 0.253NS 0.095NS 
Gracilinanus agilis 108 39 94 0.951 0.0167 −1.016NS −3.756NS 0.069NS 
Marmosa murina 80 21 49 0.908 0.013 0.151NS 1.048NS 0.104NS 
Marmosa demerarae 40 17 67 0.915 0.015 −0.952NS 0.734NS 0.088NS 
Philander opossum 54 11 0.723 0.0022 −0.743NS −1.464NS 0.0783NS 
Species n h S Hd π D FR2 
Didelphis albiventris 15 18 0.876 0.0062 −0.40 NS −1.112NS 0.1246NS 
Didelphis marsupialis 50 0.642 0.0014 −0.378NS 0.253NS 0.095NS 
Gracilinanus agilis 108 39 94 0.951 0.0167 −1.016NS −3.756NS 0.069NS 
Marmosa murina 80 21 49 0.908 0.013 0.151NS 1.048NS 0.104NS 
Marmosa demerarae 40 17 67 0.915 0.015 −0.952NS 0.734NS 0.088NS 
Philander opossum 54 11 0.723 0.0022 −0.743NS −1.464NS 0.0783NS 

NS, not significant.

Table 3.

Posterior distributions of parameter estimates from IM, including the mode number of migrants per generation (m) between populations from either Araguaia river bank, time of divergence (t) in millions of years before present (Myr BP), and the effective population size of both populations as well as of the ancestral population (qEAR, qWAR, qIslands, and qA, respectively), with 95% HPD

Species Populations ma (95% HPD) t (95% HPD) qEAR (95% HPD) qWAR (95% HPD) qIslands (95% HPD) qA (95% HPD) 
Didelphis albiventris EAR/WAR 3.11 (0.95–19.45) 0.645 (0.031–1.389) 374819 (84849–1552824) 485206 (79906–1579185) — 663142 (38717–1594013) 
Didelphis marsupialis EAR/Islands 1.29 (0.54–17.71) 0.732 (0.042–1.394) 60589 (9129–290001) — 20417 (2489–162181) 163177 (8133–323533) 
WAR/Islands 0.95 (0.41–16.95) 0.717 (0.024–1.394) — 181474 (35625–454523) 29887 (4542–146566) 240770 (13150–465999) 
EAR/WAR 1.61 (0.71–19.01) 0.702 (0.025–1.392) 38214 (9389–180152) 164430 (29479–416861) — 220331 (12010–425596) 
Philander opossum EAR/WAR 1.74 (0.615–9.115) 0.131 (0.061–1.394) 49571 (16811–138369) 487524 (117678–841853) — 407348 (18535–839266) 
Marmosa murina EAR/WAR 0.001 (0.001– 0.123) 0.451 (0.229–1.378) 228618 (107940–460589) 543723 (311752–933916) — 482043 (20783–1286565) 
Marmosa demerarae EAR/WAR 0.047 (0.007– 0.641) 0.134 (0.050–0.419) 96622 (30095–270860) 593992 (245516–1582395) — 1500028 (460937–3049159) 
Species Populations ma (95% HPD) t (95% HPD) qEAR (95% HPD) qWAR (95% HPD) qIslands (95% HPD) qA (95% HPD) 
Didelphis albiventris EAR/WAR 3.11 (0.95–19.45) 0.645 (0.031–1.389) 374819 (84849–1552824) 485206 (79906–1579185) — 663142 (38717–1594013) 
Didelphis marsupialis EAR/Islands 1.29 (0.54–17.71) 0.732 (0.042–1.394) 60589 (9129–290001) — 20417 (2489–162181) 163177 (8133–323533) 
WAR/Islands 0.95 (0.41–16.95) 0.717 (0.024–1.394) — 181474 (35625–454523) 29887 (4542–146566) 240770 (13150–465999) 
EAR/WAR 1.61 (0.71–19.01) 0.702 (0.025–1.392) 38214 (9389–180152) 164430 (29479–416861) — 220331 (12010–425596) 
Philander opossum EAR/WAR 1.74 (0.615–9.115) 0.131 (0.061–1.394) 49571 (16811–138369) 487524 (117678–841853) — 407348 (18535–839266) 
Marmosa murina EAR/WAR 0.001 (0.001– 0.123) 0.451 (0.229–1.378) 228618 (107940–460589) 543723 (311752–933916) — 482043 (20783–1286565) 
Marmosa demerarae EAR/WAR 0.047 (0.007– 0.641) 0.134 (0.050–0.419) 96622 (30095–270860) 593992 (245516–1582395) — 1500028 (460937–3049159) 

aReported as modes.

Figure 3.

Demographic scenarios of 6 didelphid marsupials from central Brazil reconstructed from mt-Cytb sequences using the program BEAST. Black solid curves indicate changes in effective population size and gray shadows indicate upper and lower 95% confidence interval.

Figure 3.

Demographic scenarios of 6 didelphid marsupials from central Brazil reconstructed from mt-Cytb sequences using the program BEAST. Black solid curves indicate changes in effective population size and gray shadows indicate upper and lower 95% confidence interval.

Didelphis marsupialis

The BI analyses based on 13 haplotypes showed 3 well-supported clades for D. marsupialis: the western Amazonia clade (WAM), the Guiana Shield clade (GS), and the Amazonia-Cerrado ecotone clade (ECO) (Figure 4). The WAM clade included all haplotypes from Juruá River, western Amazonia, with exception of a sample from locality 15 (Peru) (Figure 4). The GS clade is represented by haplotypes from the Guyana (Figure 4). The ECO clade is represented mainly by haplotypes from the Amazonia-Cerrado ecotone, especially from mid-Araguaia River, but also a sample from Tocantins River (Figure 4). The relation between these 3 clades is not resolved (Figure 4a), but their tMRCA dates to the Pleistocene at 0.44–1.65 Myr BP (Table 1, Supplementary Figure S1 online).

Figure 4.

(a) Bayesian phylogenetic reconstruction of mt-Cytb haplotypes of Didelphis marsupialis (GTR+I model used). Numbers in circle represent nodes for which time to most recent common ancestor are shown in Table 1; percentage values correspond to pairwise genetic distances between clades; and asterisks indicate BPP ≥0.95; numbers in the tip of branches are sampling localities and acronyms are haplotype designations; clades are designated by ECO (Amazonia-Cerrado ecotone), GS (Guiana Shield), and WAM (western Amazonia). (b) MJ network of D. marsupialis of the ECO clade, central Brazil. Length of connecting branches corresponds to the nucleotide substitutions; size of the circles is proportional to the number of individuals sharing each haplotype; white corresponds to the eastern riverbank of the Araguaia River, black to the western riverbank, and gray to the islands of the Araguaia River. (c) Map of sampling localities of D. marsupialis, with the representation of the clades.

Figure 4.

(a) Bayesian phylogenetic reconstruction of mt-Cytb haplotypes of Didelphis marsupialis (GTR+I model used). Numbers in circle represent nodes for which time to most recent common ancestor are shown in Table 1; percentage values correspond to pairwise genetic distances between clades; and asterisks indicate BPP ≥0.95; numbers in the tip of branches are sampling localities and acronyms are haplotype designations; clades are designated by ECO (Amazonia-Cerrado ecotone), GS (Guiana Shield), and WAM (western Amazonia). (b) MJ network of D. marsupialis of the ECO clade, central Brazil. Length of connecting branches corresponds to the nucleotide substitutions; size of the circles is proportional to the number of individuals sharing each haplotype; white corresponds to the eastern riverbank of the Araguaia River, black to the western riverbank, and gray to the islands of the Araguaia River. (c) Map of sampling localities of D. marsupialis, with the representation of the clades.

The MJ network was constructed with 50 samples of the ECO clade (Table 2). Haplotype HDm1 is the ancestor of the Araguaia River samples and it is shared by specimens from both riverbanks (Figure 4b). Haplotype HDm3 is also shared by specimens from both riverbanks, along with specimens from the islands of the Araguaia River (Figure 4b). Haplotype HDm2 was only sampled in islands, HDm4 in the eastern bank, and HDm5 in the western bank of the Araguaia River (Figure 4b). The IM results of the ECO clade suggested lower gene flow than D. albiventris, but still substantial ranging from 0.95 (95% HPD = 0.41–16.95) for western river bank and islands to 1.61 (95% HPD = 0.71–19.01) for eastern and western river banks (Table 3, Supplementary Figure S2 online). Although the divergence time t did not converged well for the 3 pairwise combinations, the best estimate of divergence time between the eastern and western river banks was 0.702 Myr BP (95% HPD = 0.025–1.392 Myr BP, Table 3). The eastern population (qEAR) had a lower inferred size than the ancestral (qA) and the western (qWAR) populations (Table 3, Supplementary Figure S3 online). Neutrality tests were not significant (Table 2), that is, they did not reject the constant population model. BSP for the ECO clade also support a constant population (Figure 3).

Philander opossum

The BI analyses based on 15 haplotypes showed 2 well-supported clades for P. opossum: the western Amazonia clade (WAM) and the Araguaia River clade (AR) (Figure 5a). The WAM clade included all haplotypes from Juruá River, western Amazonia (Figure 5). The AR clade is represented by haplotypes from the Araguaia River, including Amazonia-Cerrado ecotone and Cerrado domains (Figure 5). Haplotypes from the Pantanal (PAN) formed a polytomy with those from the Araguaia River (AR), central Brazil clade (Figure 5).

Figure 5.

(a) Bayesian phylogenetic reconstruction of mt-Cytb haplotypes of Philander opossum (GTR+G model used). Numbers in circle represent nodes for which time to most recent common ancestor are shown in Table 1; percentage values correspond to pairwise genetic distances between clades; and asterisks indicate BPP ≥0.95; numbers in the tip of branches are sampling localities and acronyms are haplotype designations; clades are designated by WAM (western Amazonia), AR (Araguaia River), and PAN (Pantanal). (b) MJ network of P. opossum of the AR clade, central Brazil. Length of connecting branches corresponds to the nucleotide substitutions; size of the circles is proportional to the number of individuals sharing each haplotype; white corresponds to the eastern riverbank of the Araguaia River, and black to the western riverbank. (c) Map of sampling localities of P. opossum, with the representation of the clades.

Figure 5.

(a) Bayesian phylogenetic reconstruction of mt-Cytb haplotypes of Philander opossum (GTR+G model used). Numbers in circle represent nodes for which time to most recent common ancestor are shown in Table 1; percentage values correspond to pairwise genetic distances between clades; and asterisks indicate BPP ≥0.95; numbers in the tip of branches are sampling localities and acronyms are haplotype designations; clades are designated by WAM (western Amazonia), AR (Araguaia River), and PAN (Pantanal). (b) MJ network of P. opossum of the AR clade, central Brazil. Length of connecting branches corresponds to the nucleotide substitutions; size of the circles is proportional to the number of individuals sharing each haplotype; white corresponds to the eastern riverbank of the Araguaia River, and black to the western riverbank. (c) Map of sampling localities of P. opossum, with the representation of the clades.

The MJ network was constructed using 54 samples of the AR clade (Figure 5b, Table 2). The haplotype HPo1 is the ancestor of the mid-Araguaia River samples and it is shared by specimens from both riverbanks (Figure 5b). Haplotypes HPo6 and HPo7 are also shared by specimens from both riverbanks (Figure 5b). The IM results of the AR clade suggested some gene flow between eastern and western river banks, with the posterior distributions of migration rate m peaking at 1.74 (95% HPD = 0.615–9.115, Table 3). The time of divergence t did not converged well (see Supplementary Figure S2 online), probably due to the high values of gene flow which could mask a signal of population separation. The eastern population (qEAR) had a lower inferred size than the ancestral (qA) and the western (qWAR) populations (Table 3, Supplementary Figure S3 online). Neutrality tests were not significant (Table 2), that is, they did not reject the constant population model. BSP for the AR clade also support a constant population (Figure 3).

Gracilinanus agilis

The BI analyses based on 65 haplotypes showed 6 well-supported clades for G. agilis: the east clade (EAC), the north-eastern clade (NEC), the western clade (WCE), the southern clade (SC), the central Brazil plateaux clade (CBP), and the central-east clade (CEC, Figure 6a). The east clade (EAC) included almost all haplotypes from the moist forest enclaves in the east of the São Francisco River, with exception of 2 haplotypes from locality 65 (Rio de Contas, state of Bahia, Figure 6). The north-eastern clade (NEC) is represented by haplotypes from the north-eastern Caatinga and Cerrado, which are located in the west margin of São Francisco River (Figure 6). The western clade (WCE) included haplotypes from Amazonia-Cerrado ecotone and western Cerrado (Figure 6). The southern clade (SC) included haplotypes from southern Cerrado and Pantanal (Figure 6). The central Brazil plateaux clade (CBP) and central-east clade (CEC) included haplotypes from central Brazil plateaux and from central-east Cerrado, respectively (Figure 6)

Figure 6.

(a) Bayesian phylogenetic reconstruction of mt-Cytb haplotypes of Gracilinanus agilis (GTR+I+G model used). Numbers in circle represent nodes for which time to most recent common ancestor are shown in Table 1; percentage values correspond to pairwise genetic distances between clades; and asterisks indicate BPP ≥0.95; numbers in the tip of branches are sampling localities and acronyms are haplotype designations; clades are designated by the SC (southern clade), CBP (central Brazil Plateau clade), CEC (central-east clade), WCE (western clade), NEC (north-eastern clade), and EAC (east clade). (b) MJ network of G. agilis of WCE clade, central Brazil. Length of connecting branches corresponds to the nucleotide substitutions; size of the circles is proportional to the number of individuals sharing each haplotype; white corresponds to the eastern riverbank of the Araguaia River, and black to the western riverbank. (c) Map of sampling localities of G. agilis, with the representation of the clades.

Figure 6.

(a) Bayesian phylogenetic reconstruction of mt-Cytb haplotypes of Gracilinanus agilis (GTR+I+G model used). Numbers in circle represent nodes for which time to most recent common ancestor are shown in Table 1; percentage values correspond to pairwise genetic distances between clades; and asterisks indicate BPP ≥0.95; numbers in the tip of branches are sampling localities and acronyms are haplotype designations; clades are designated by the SC (southern clade), CBP (central Brazil Plateau clade), CEC (central-east clade), WCE (western clade), NEC (north-eastern clade), and EAC (east clade). (b) MJ network of G. agilis of WCE clade, central Brazil. Length of connecting branches corresponds to the nucleotide substitutions; size of the circles is proportional to the number of individuals sharing each haplotype; white corresponds to the eastern riverbank of the Araguaia River, and black to the western riverbank. (c) Map of sampling localities of G. agilis, with the representation of the clades.

The relation between most of these clades is not resolved, but their tMRCA dates to the Pliocene/Pleistocene at 1.15–3.32 Myr BP (node 3, Figure 6a, Table 1). The average genetic distance between clades varied from 4.7±0.6% (EAC and the rest of the clades) to 2.5±0.4% (WCE and CBP). The tMRCA of central Brazilian clades (WCE + CBP + CEC + SC) was estimated during the Pleistocene at 0.71–1.94 Myr BP (Table 1, Supplementary Figure S1 online).

The MJ network was constructed using 56 samples of the WCE clade (Figure 6b). The haplotype HGa4 is centrally located relatively to samples from the eastern bank of the Araguaia River (Figure 6b). The haplotype HGa11 was only sampled in the western bank of the Araguaia River (Figure 6b). Neutrality tests were not significant (Table 2), that is, they did not reject the constant population model. BSP for the WCE clade also support a constant population (Figure 3).

Marmosa murina

The BI analyses based on 36 haplotypes showed 5 well-supported clades for M. murina: the Atlantic Forest clade (AF), the Tapajós River clade (TR), the Araguaia-Tocantins interfluve clade (ATI), the western Araguaia River clade (WAR) and the eastern Araguaia River clade (EAR) (Figure 7a). Haplotypes from the Guiana Shield (GS) and northern Amazonia (NAM) form a polytomy with the remaining haplotypes (Figure 7a).

Figure 7.

(a) Bayesian phylogenetic reconstruction of mt-Cytb haplotypes of Marmosa murina (HKY+G model used). Numbers in circle represent nodes for which time to most recent common ancestor are shown in Table 1; percentage values correspond to pairwise genetic distances between clades; and asterisks indicate BPP ≥0.95; numbers in the tip of branches are sampling localities and acronyms are haplotype designations; clades are designated by EAR (eastern Araguaia River), WAR (western Araguaia River), ATI (Araguaia-Tocantins interfluve), TR (Tapajos River), AF (Atlantic Forest), NAM (northern Amazonia), and GS (Guiana Shield). (b) MJ network of M. murina of EAR, WAR, and ATI clades, central Brazil. Length of connecting branches corresponds to the nucleotide substitutions; size of the circles is proportional to the number of individuals sharing each haplotype; white corresponds to the eastern riverbank of the Araguaia River, and black to the western riverbank. (c) Map of sampling localities of M. murina, with the representation of the clades.

Figure 7.

(a) Bayesian phylogenetic reconstruction of mt-Cytb haplotypes of Marmosa murina (HKY+G model used). Numbers in circle represent nodes for which time to most recent common ancestor are shown in Table 1; percentage values correspond to pairwise genetic distances between clades; and asterisks indicate BPP ≥0.95; numbers in the tip of branches are sampling localities and acronyms are haplotype designations; clades are designated by EAR (eastern Araguaia River), WAR (western Araguaia River), ATI (Araguaia-Tocantins interfluve), TR (Tapajos River), AF (Atlantic Forest), NAM (northern Amazonia), and GS (Guiana Shield). (b) MJ network of M. murina of EAR, WAR, and ATI clades, central Brazil. Length of connecting branches corresponds to the nucleotide substitutions; size of the circles is proportional to the number of individuals sharing each haplotype; white corresponds to the eastern riverbank of the Araguaia River, and black to the western riverbank. (c) Map of sampling localities of M. murina, with the representation of the clades.

The monophyly of each clade from both river banks of the Araguaia River (WAR and EAR) and from the ATI is supported with high BPP (1.00, Figure 7a). The tMRCA of these clades dates to the Pleistocene at 0.39–1.32 Myr BP (Table 1), and the average genetic distance between them is 1.9±0.4%.

The MJ network was constructed using 80 samples of WAR, EAR, and ATI clades (Figure 7b, Table 2). Our results also suggested that ancestors of these clades have not been sampled yet and haplotypes are exclusive from each river bank (WAR and EAR) (Figure 7b). Haplotypes from ATI are divergent and separated by 8 mutations, resulting in a MJ network with large branches, and with no intermediary haplotypes (Figure 7b). The IM results of WAR and EAR clades suggested very low gene flow between eastern and western river banks, with the posterior distributions of migration rate m peaking at 0.001 (95% HPD = 0.001–0.123, Table 3, Supplementary Figure S2 online). The posterior distribution of time of divergence t peaked at 6.31 (95% HPD = 3.21–19.01, Supplementary Figure S2 online), which converted to time in years suggests that the populations split occurred about 0.451 Myr BP (95% HPD = 0.229–1.378 Myr BP Table 3). The eastern population (qEAR) had a lower inferred size than the ancestral (qA) and the western (qWAR) populations (Table 3, Supplementary Figure S3 online). Neutrality tests were not significant (Table 2), that is, they did not reject the constant population model. BSP for the EAR and WAR clades shows a slightly population decline to the present (Figure 3).

Marmosa demerarae

The BI analyses based on 37 haplotypes showed 6 well-supported clades for M. demerarae: the western Amazonia clade (WAM), the northern Amazonia clade (NAM), the Guiana Shield clade (GS), the Amazonia-Cerrado ecotone clade (ECO), the eastern Tocantins River clade (EToR), and the Atlantic Forest clade (AF) (Figure 8a). The WAM clade is a divergent clade (5.7±0.7%), including samples from southern Amazon River to the Teles Pires River (Figure 8c). The AF clade is represented by almost all haplotypes from the Atlantic Forest domain, with exception of a haplotype from locality 60 (São José de Lages, state of Alagoas) in the northern Atlantic Forest (Figure 8). This haplotype is the only one located west of the São Francisco River, while remaining samples from the Atlantic Forest are located east of this river (Figure 8c). The ECO clade includes haplotypes from the Amazonia-Cerrado ecotone, mainly from both Araguaia riverbanks but also a haplotypes from the Xingu River (Figure 8c). The EToR clade includes haplotypes from the eastern bank of the lower Tocantins River in the Amazonia-Cerrado ecotone (HMd15 and HMd16) and in the Cerrado (HMd17). Haplotype HMd12 from the western bank of the lower Tocantins River (locality 30) is very distant (4.2±0.7%) from the rest of the samples of central Brazil, and formed a polytomy with AF, ECO, and EToR clades (Figure 8).

Figure 8.

(a) Bayesian phylogenetic reconstruction of mt-Cytb haplotypes of Marmosa demerarae (HKY+I+G model used). Numbers in circle represent nodes for which time to most recent common ancestor are shown in Table 1; percentage values correspond to pairwise genetic distances between clades; and asterisks indicate BPP ≥0.95; numbers in the tip of branches are sampling localities and acronyms are haplotype designations; clades are designated by WAM (western Amazonia), NAM (northern Amazonia), GS (Guiana Shield), ECO (Amazonia-Cerrado ecotone), EToR (eastern Tocantins River), and AF (Atlantic Forest). (b) MJ network of M. demerarae of the ECO (including samples from eastern Araguaia River, EAR) and EToR clades, central Brazil. Length of connecting branches corresponds to the nucleotide substitutions; size of the circles is proportional to the number of individuals sharing each haplotype; white corresponds to the eastern riverbank of the Araguaia River, and black to the western riverbank. (c) Map of sampling localities of M. demerarae, with the representation of the clades.

Figure 8.

(a) Bayesian phylogenetic reconstruction of mt-Cytb haplotypes of Marmosa demerarae (HKY+I+G model used). Numbers in circle represent nodes for which time to most recent common ancestor are shown in Table 1; percentage values correspond to pairwise genetic distances between clades; and asterisks indicate BPP ≥0.95; numbers in the tip of branches are sampling localities and acronyms are haplotype designations; clades are designated by WAM (western Amazonia), NAM (northern Amazonia), GS (Guiana Shield), ECO (Amazonia-Cerrado ecotone), EToR (eastern Tocantins River), and AF (Atlantic Forest). (b) MJ network of M. demerarae of the ECO (including samples from eastern Araguaia River, EAR) and EToR clades, central Brazil. Length of connecting branches corresponds to the nucleotide substitutions; size of the circles is proportional to the number of individuals sharing each haplotype; white corresponds to the eastern riverbank of the Araguaia River, and black to the western riverbank. (c) Map of sampling localities of M. demerarae, with the representation of the clades.

The MJ network was constructed using 40 samples of ECO and EToR clades (Figure 8b, Table 2). Haplotypes HMd1, HMd3–HMd5, and HMd7 from the eastern bank of the mid-Araguaia River (EAR), have as ancestor the haplotype HMd14 from the western bank of the lower Araguaia River (locality 28, Figure 8bc). The EToR clade and the haplotype HMd12 are separated by 25 and 28 mutations, respectively, from the rest of the samples, which result in a MJ network with very large branches, and with no intermediary haplotypes (Figure 8b). The IM results suggested very limited gene flow between eastern and western river banks, with the posterior distributions of migration rate m peaking at 0.047 (95% HPD = 0.007–0.641, Table 3, Supplementary Figure S2 online). Although the divergence time t did not converged well, there was a peak at 1.88 (95% HPD = 0.693–5.859, Supplementary Figure S2 online), which converted to time in years suggests that the populations started diverging about 0.134 Myr BP (95% HPD = 0.050–0.419 Myr BP, Table 3). The eastern population (qEAR) had a lower inferred size than the ancestral (qA) and the western (qWAR) populations (Table 3, Supplementary Figure S3 online). Neutrality tests were not significant (Table 2), that is, they did not reject the constant population model. BSP for the ECO clade also support a constant population (Figure 3).

Discussion

Comparative Phylogeography and Species Ecology

The 6 didelphid marsupials analyzed in this work revealed different patterns of differentiation across the mid-Araguaia River. The 3 larger species, D. albiventris (body size, i.e., head and body length: 265–363mm), D. marsupialis (310–460mm) and P. opossum (180–260mm), revealed extensive intraspecific haplotype sharing across the river (Figures 2b, 4b, and 5b). Consistent with this haplotype sharing, the IM results also confirmed gene flow among Araguaia river banks for the 2 latter species (Table 3).

In contrast, the 3 smaller marsupial species (G. agilis, M. murina, and M. demerarae, body sizes: 75–108, 98–155, and 110–190mm, respectively) showed no haplotype sharing across the mid-Araguaia River. For G. agilis we did not observe haplotype sharing across the river (Figure 6b), but the geographic gap between our samples in this region makes the role of the river as a barrier less obvious. Therefore, only increasing sampling effort in the western river bank will help to proper access the role of the Araguaia River in the differentiation of this species. For M. demerarae, although the paraphyletic relationship was not supported in the Bayesian analysis (Figure 8a), the observed ancestral relationships in the MJ network (Figure 8b) revealed that the western bank population was most likely the origin from which the eastern bank was colonized. This pattern is also supported by the observation that for M. demerarae, as well as for D. marsupialis, P. opossum, and M. murina, the eastern population had lower inferred size than the ancestral and the western populations. Although tentative this pattern could be explained by the fact that these 4 species are either typically Amazonian species (D. marsupialis and P. opossum) or that occur mainly in rainforests (M. murina and M. demerarae), and the western bank of the mid-Araguaia might represent a previous distribution limit of these species.

Haplotypes of M. murina on opposite riverbanks formed monophyletic clades that are genetically distant from each other (Figure 7a). Although reciprocal monophyly, which would confirm the formation of a complete barrier to an existing taxon (Patton et al. 2000), was not observed, we hypothesized that this barrier could have been formed by steps. Indeed, the Araguaia River is acting as barrier to gene flow of populations on opposite river banks as indicated by IM results of M. murina (and also M. demeraraeTable 3). Additionally, IM results also suggest that eastern and western populations of Araguaia River most likely split about 0.451 Myr BP (0.229–1.379 Myr BP) for M. murina, and about 0.134 Myr BP (95% HPD = 0.050–0.419 Myr BP) for M. demerarae, during the Pleistocene. This divergence estimates encompass the age of Araguaia River sediments formation (roughly 0.240 Myr BP, Valente and Latrubesse 2012). Therefore, we hypothesize that the Araguaia might have formed a barrier to populations of M. murina in the mid-Araguaia (WAR and EAR) and later to the populations in the lower Araguaia (ATI), but the 2 populations in the eastern Araguaia River (EAR and ATI) have also limited gene flow among them (see Other phylogeographic patterns to further discussion on ATI divergence).

Moreover, it should be noted that the Araguaia River marks the transition between 2 different biomes: the Amazonia and the Cerrado. Therefore, it is possible that the differentiation observed across the river and in this region could be the result of local adaptation to 2 different biomes, with reduced effective gene flow due to reduced fitness of admixed individuals (e.g., Almeida et al. 2007; Miranda et al. 2007, 2009). If genetic differentiation was caused by adaptation to different biomes we would expect sister clades related with adjacent but distinct habitats (Moritz et al. 2000; Smith et al. 2005). However, limited difference in habitat between the 2 riverbanks allows us to discard this hypothesis in explaining the intraspecific genetic structure of M. murina and M. demerarae. Additionally, a signature of demographic stability detected for didelphid species along the Araguaia River do not support the diversification in refuges. Although the importance of Pleistocene climatic oscillations as a driver of population divergence in central Brazil is still in debate, the scenario of limited or no demographic expansion has also been reported to small mammals in the western Amazonia, contradicting the refugia model of diversification (Matocq et al. 2000; Lessa et al. 2003). On the other hand, rivers of the central Brazilian plateux have been shown to be stronger historical barriers in promoting differentiation (e.g., Maldonado-Coelho et al. 2013), as corroborated by our results.

Differences in phylogeography among species are likely to be linked to differences in the life history (e.g., mating system and reproduction success) and ecological traits (e.g., locomotion mode and habitat preferences) (Matocq et al. 2000; Rocha et al. 2014). Although body size separates the 2 groups of didelphid species, this morphological feature is not likely to be the main cause of the observed genetic differences in the mid-Araguaia River. Small rodents (body size ranging from 90 to 125mm) showed distinct patterns of genetic differentiation in this river (Rocha et al. 2011, 2014). In addition to body size, habitat preferences and overall distribution differ among didelphid species (Patton et al. 2000; Gardner 2008; Rocha et al. 2011). These ecological differences probably have a great impact on their genetic structure. Marmosa murina is known to be associated with upland forests and avoids gallery forests that flood (Patton et al. 2000; Ramos Pereira et al. 2013). Marmosa demerarae and G. agilis are not upland forest specialists, since they also use seasonally flooded forests to some extent (Rocha et al. 2011; Ramos Pereira et al. 2013). These arboreal murine mouse opossums search the understory and forest floor for invertebrate prey (Emmons and Feer 1997; Gardner 2008); a foraging strategy impossible to accomplish during the flooded season, which may explain their preference for dry areas (Ramos Pereira et al. 2013). Contrastingly, all of D. albiventris, D. marsupialis, and P. opossum are well distributed both in upland and seasonally flooded forests (Rocha et al. 2011), which have probably enhanced cross-river gene flow. Indeed, previous authors have pointed out that upland forest specialists are more prone to exhibit riverine diversification than specialists of floodplain forests (Moritz et al. 2000; Patton et al. 2000; Rocha et al. 2014), being the latter capable of colonizing river islands and crossing rivers due life-history attributes such as high dispersal rates (Aleixo 2006).

The Importance of River Geomorphology in Biogeography

River basin geomorphology may be an important factor determining the importance of the river as a phylogeographic barrier. Meandering rivers, where the location of the river bed changes over time, as the Juruá River (Latrubesse 2008), may facilitate the exchange of individuals by repeatedly bringing vicariant populations into contact (Moritz et al. 2000; Patton et al. 2000). Contrastingly, anabranching rivers with low sinuosities, as the rivers in the Brazilian Shield, including Araguaia River (Latrubesse 2008), may be stronger barriers simply because they are more permanent physical obstacles. These predictions were corroborated by previous studies at the Juruá River, which was not a barrier for most of the sampled vertebrates (e.g., Gascon et al. 2000; Matocq et al. 2000; Patton et al. 2000), and at Tocantins and Teles Pires Rivers, which were found to be historical barriers in the differentiation of several bird species (Bates et al. 2004; Maldonado-Coelho et al. 2013). They are also confirmed by our results at the mid-Araguaia River, although we find the phylogeographical impact of the river to be species-dependent (see also Rocha et al. 2011, 2014; Faria et al. 2013a). Our comparative results suggest that semi-permanent alluvial islands characteristic from anabranching rivers may be important stepping-stones facilitating the gene flow of large didelphids, as well as small rodents (Rocha et al. 2014), which inhabit seasonally flooded forests. Indeed, haplotype sharing and gene flow among populations of D. marsupialis of both Araguaia river banks and of the islands was confirmed by our data.

Other Phylogeographic Patterns

Phylogeographic analysis of M. demerarae suggests that the dispersion of this species in central Brazil occurred throughout the gallery forests of the Araguaia and Tocantins River. Specifically, the observed genetic pattern revealed that specimens from the eastern bank of the Tocantins River are genetically similar despite the large geographic distance, while most of the specimens from the eastern and western bank of the Araguaia River belong to the same cluster (Figure 8c). This finding supports previous studies that have shown that these gallery forest act as dispersal corridors for rainforest mammal species throughout the Cerrado (e.g., Johnson et al. 1999; Costa 2003).

Very divergent haplotypes of M. demerarae and M. murina were recorded in the ATI. Previously, a potentially new species of the subgenus Marmosa (Micoureus) was recorded circa 300 km east of this interfluve (Costa 2003), revealing the importance of this region to the diversification within this subgenus. In addition to the river interfluve, palynological records of the Serra dos Carajás, state of Pará, showed that this region went through several alternating periods dominated by arboreal and herbaceous savanna vegetations (Iriondo and Latrubesse 1994; Behling et al. 2010). However previous authors fail in support bird differentiation in glacial forest refuges, instead they found evidences of riverine barrier caused by the Tocantins River (Maldonado-Coelho et al. 2013). Further intensive sampling in this region will help to uncover the ATI contribution to the diversification of these didelphid marsupials.

Finally, our phylogeographic analyses of G. agilis revealed highly supported clades within Cerrado of central Brazil. Faria et al. (2013b) discussed the importance of the Serra Geral do Goiás to the diversification of this species. Although the geological origin of these mountain chains dates to the Cretaceous (Villela and Nogueira 2011), before the diversification of this genus, the western clades (SC, CBP, CEC, and WCE) and northern-east clade (NEC) are separated by the Serra Geral do Goiás (Faria et al. 2013b, and also shown in the present study). Moreover, 7 geographical groups within the Cerrado were recognized based on its floristic composition (Ratter et al. 2003). The western clades of G. agilis are not completely congruent with the geographical limits of the floristic groups (Ratter et al. 2003). However, vegetation pattern reflects changes that occurred during the Tertiary and Quaternary periods (Ratter et al. 2003), and the extant diversity of G. agilis in Cerrado dates at 1.26 Myr BP (0.71–1.94 Myr BP). Thus, suggesting that vegetation fluctuations occurred during this period may have influenced the intraspecific genetic structure of this species (see also Nascimento et al. 2013).

Supplementary Material

Supplementary material can be found at http://www.jhered.oxfordjournals.org/

Funding

PhD grants by Portuguese Science Foundation (FCT—Program POPH-QREN, refs: SFRH/BD/24767/2005 and SFRH/BD/23191/2005 to R.G.R. and E.F.); Research grants from Fundação de Amparo à Pesquisa e Inovação do Espírito Santo (FAPES, Brazil to R.G.R. and L.P.C.). PhD scholarship from Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES—PhD scholarship to A.C.L.); Fellowship from Conselho Nacional de Desenvolvimento Científico e Tenológico (CNPq, Brazil to L.P.C.); and European Funds through COMPETE, by the FCT within project PEst-C/MAR/LA0017/2013 and by CNPq (Rede BioM.A., Processo: 457491/2012–4).

Acknowledgments

We thank Rafael Leite and Jim Patton who kindly shared samples and sequences of didelphid marsupials; Juliana Justino who provided crucial help in the laboratory; Yuri L.R. Leite who provided helpful comments throughout the work. We also thank anonymous reviewers for helpful comments and suggestions that improved this manuscript. Universidade Federal do Tocantins/ Fundação de Apoio Científico e Tecnológico do Tocantins (UFT/FAPTO) and Ecotropical (a partnership between Instituto Ecológica and Universidade de Aveiro) provided logistic support. Samples were collected under the permits of Federal (Instituto Chico Mendes de Conservação da Biodiversidade-CMBio, permits 200/2006 and 14307-1) and State (NATURATINS, permits: 019/2006 and 001/2008) conservancy agencies.

References

Aleixo
A
.
2006
.
Historical diversification of a terra-firme forest bird superspecies: a phylogeographic perspective on the role of different hypotheses of Amazonian diversification
.
Evolution
 .
58
:
1303
1317
.
Almeida
FC
Bonvicino
CR
Cordeiro-Estrela
P
.
2007
.
Phylogeny and temporal diversification of Calomys (Rodentia, Sigmodontinae): implications for the biogeography of an endemic genus of the open/dry biomes of South America
.
Mol Phylogenet Evol
 .
42
:
449
466
.
Antonelli
A
Quijada–Mascarenas
A
Crawford
AJ
Bates
JM
Velazco
PM
Wuster
W
.
2010
.
Molecular studies and phylogeography of Amazonian tetrapods and their relation to geological and climatic models
. In:
Hoorn
C
Wesselingh
F
, editors.
Amazonia, landscape and species evolution: a look into the past
 .
West Sussex (UK)
:
Wiley–Blackwell
. p.
386
404
.
Arbogast
BS
Kenagy
GJ
.
2001
.
Comparative phylogeography as an integrative approach to historical biogeography
.
J Biogeogr
 .
28
:
819
825
.
Ayres
JM
Clutton-Brock
TH
.
1992
.
River boundaries and species range size in Amazonian primates
.
Am Nat
 .
140
:
531
537
.
Ballard
JW
Whitlock
MC
.
2004
.
The incomplete natural history of mitochondria
.
Mol Ecol
 .
13
:
729
744
.
Bandelt
HJ
Forster
P
Röhl
A
.
1999
.
Median-joining networks for inferring intraspecific phylogenies
.
Mol Biol Evol
 .
16
:
37
48
.
Baker CS.
2013
.
Journal of Heredity adopts joint data archiving policy
.
J Hered
 .
104
:
1
.
Bates
JM
Haffer
J
Grismer
E
.
2004
.
Avian mitochondrial DNA sequence divergence across a headwater stream of the Rio Tapajós, a major Amazonian river
.
J Ornithol
 .
145
:
199
205
.
Behling
H
Bush
M
Hooghiemstra
H
.
2010
.
Biotic development of Quaternary Amazonia: a palynological perspective
. In:
Hoorn
C
Wesselingh
F
, editors.
Amazonia, landscape and species evolution: a look into the past
 .
West Sussex (UK)
:
Wiley–Blackwell
. p.
335
345
.
Bezerra
AMR
Carmignotto
AP
Rodrigues
FHG
.
2009
.
Small non-volant mammals of an ecotone region between the Cerrado hotspot and the Amazonian rainforest, with comments on their taxonomy and distribution
.
Zool Stud
 .
48
:
861
874
.
Bouckaert
R
Heled
J
Kühnert
D
Vaughan
T
Wu
CH
Xie
D
Suchard
MA
Rambaut
A
Drummond
AJ
.
2014
.
BEAST 2: a software platform for Bayesian evolutionary analysis
.
PLoS Comput Biol
 .
10
:
e1003537
.
Bruford
MW
Hanotte
O
Brookfield
JFY
Burke
T
.
1992
.
Single-locus and DNA fingerprinting
. In:
Hoelzel
AR
, editor.
Molecular genetic analyses of populations. A pratical approach
 .
Oxford
:
IRL Press
. p.
225
269
.
Cohen
KM
Finney
SC
Gibbard
PL
Fan
JX
.
2013
.
The ICS International Chronostratigraphic Chart
.
Episodes
 .
36
:
199
204
.
Costa
LP
.
2003
.
The historical bridge between the Amazon and the Atlantic Forest of Brazil: a study of molecular phylogeography with small mammals
.
J Biogeogr
 .
30
:
71
86
.
Drummond
AJ
Ho
SY
Phillips
MJ
Rambaut
A
.
2006
.
Relaxed phylogenetics and dating with confidence
.
PLoS Biol
 .
4
:
699
710
.
Emmons
LH
Feer
F
.
1997
.
Neotropical rainforest mammals, second edition
 .
Chicago (IL)
:
University of Chicago Press
.
Faria
MB
de Oliveira
JA
Bonvicino
CR
.
2013a
.
Filogeografia de populações brasileiras de Marmosa (Marmosa) murina (Didelphimorphia, Didelphidae)
.
Rev Bras Biol
 .
21
:
27
52
.
Faria
MB
Nascimento
FF
Oliveira
JA
Bonvicino
CR
.
2013b
.
Biogeographic determinants of genetic diversification in the mouse opossum Gracilinanus agilis (Didelphimorphia: Didelphidae)
.
J Hered
 .
104
:
613
626
.
Fu
YX
.
1997
.
Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection
.
Genetics
 .
147
:
915
925
.
Gardner
AL
.
2008
.
Mammals of South America, vol. 1. Marsupials, xenarthrans, shrews, and bats
 .
Chicago (IL)
:
Chicago University Press
.
Gascon
C
Malcolm
JR
Patton
JL
da Silva
MN
Bogart
JP
Lougheed
SC
Peres
CA
Neckel
S
Boag
PT
.
2000
.
Riverine barriers and the geographic distribution of Amazonian species
.
Proc Natl Acad Sci U S A
 .
97
:
13672
13677
.
Hey
J
Nielsen
R
.
2004
.
Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. persimilis
.
Genetics
 .
167
:
747
760
.
Horovitz
I
Martin
T
Bloch
J
Ladevèze
S
Kurz
C
Sánchez-Villagra
MR
.
2009
.
Cranial anatomy of the earliest marsupials and the origin of opossums
.
PLoS One
 .
4
:
e8278
.
Iriondo
M
Latrubesse
EM
.
1994
.
A probable scenario for a dry climate in central Amazonia during the late Quaternary
.
Quaternary Int
 .
21
:
121
128
.
Jansa
SA
Barker
FK
Voss
RS
.
2014
.
The early diversification history of didelphid marsupials: a window into South America’s “Splendid Isolation”
.
Evolution
 .
68
:
684
695
.
Johnson
MA
Saraiva
PM
Coelho
D
.
1999
.
The role of gallery forest in the distribution of Cerrado mammals
.
Rev Bras Biol
 .
59
:
421
427
.
Kimura
M
.
1980
.
A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences
.
J Mol Evol
 .
16
:
111
120
.
Lacher
TE
Alho
CJR
.
2001
.
Terrestrial small mammal richness and habitat associations in an Amazon forest-Cerrado contact zone
.
Biotropica
 .
33
:
171
181
.
Latrubesse
EM
.
2008
.
Patterns of anabranching channels: the ultimate end-member adjustment of mega rivers
.
Geomorphology
 .
101
:
130
145
.
Latrubesse
EM
Stevaux
JC
.
2002
.
Geomorphology and environmental aspects of the Araguaia fluvial basin, Brazil
.
Z Geomorphol
 .
129
:
109
127
.
Leite
R
Rogers
D
.
2013
.
Revisiting Amazonian phylogeography: insights into diversification hypotheses and novel perspectives
.
Org Divers Evol
 .
13
:
639
664
.
Lessa
EP
Cook
JA
Patton
JL
.
2003
.
Genetic footprints of demographic expansion in North America, but not Amazonia, during the Late Quaternary
.
Proc Natl Acad Sci U S A
 .
100
:
10331
10334
.
Librado
P
Rozas
J
.
2009
.
DnaSP v5: a software for comprehensive analysis of DNA polymorphism data
.
Bioinformatics
 .
25
:
1451
1452
.
Maldonado-Coelho
M
Blake
JG
Silveira
LF
Batalha-Filho
H
Ricklefs
RE
.
2013
.
Rivers, refuges and population divergence of fire-eye antbirds (Pyriglena) in the Amazon Basin
.
J Evol Biol
 .
26
:
1090
1107
.
Matocq
MD
Patton
JL
da Silva
MN
.
2000
.
Population genetic structure of two ecologically distinct Amazonian spiny rats: separating history and current ecology
.
Evolution
 .
54
:
1423
1432
.
Meredith
RW
Janečka
JE
Gatesy
J
Ryder
OA
Fisher
CA
Teeling
EC
Goodbla
A
Eizirik
E
Simão
TL
Stadler
T
et al.  
2011
.
Impacts of the cretaceous terrestrial revolution and KPg extinction on mammal diversification
.
Science
 .
334
:
521
524
.
Miranda
GB
Andrades-Miranda
J
Oliveira
LF
Langguth
A
Mattevi
MS
.
2007
.
Geographic patterns of genetic variation and conservation consequences in three South American rodents
.
Biochem Genet
 .
45
:
839
856
.
Miranda
GB
Oliveira
LF
Andrades-Miranda
J
Langguth
A
Callegari-Jacques
SM
Mattevi
MS
.
2009
.
Phylogenetic and phylogeographic patterns in sigmodontine rodents of the genus oligoryzomys
.
J Hered
 .
100
:
309
321
.
Moritz
C
Patton
JL
Schneider
CJ
Smith
TB
.
2000
.
Diversification of rainforest faunas: an integrated molecular approach
.
Annu Rev Ecol Syst
 .
31
:
533
563
.
Nascimento
FF
Lazar
A
Menezes
AN
Durans
Ada M
Moreira
JC
Salazar-Bravo
J
D’Andrea
PS
Bonvicino
CR
.
2013
.
The role of historical barriers in the diversification processes in open vegetation formations during the Miocene/Pliocene using an ancient rodent lineage as a model
.
PLoS One
 .
8
:
e61924
.
Nicolas
V
Missoup
AD
Denys
C
Peterhans
JK
Katuala
P
Couloux
A
Colyn
M
.
2011
.
The roles of the rivers and Pleistocene refugia in shaping genetic diversity in Praomys misonnei in tropical Africa
.
J Biogeogr
 .
38
:
191
207
.
Nylander
JAA
.
2004
.
MrModeltest, version 2. Program distributed by the author
 .
Uppsala (Sweden)
:
Evolutionary Biology Centre, Uppsala University
.
Patton
JL
da Silva
MNF
Malcolm
JR
.
2000
.
Mammals of the Rio Juruá and the evolutionary and ecological diversification of Amazonia
.
Bull Am Mus Nat Hist
 .
244
:
1
306
.
Peters
JL
Gretes
W
Omland
KE
.
2005
.
Late Pleistocene divergence between eastern and western populations of wood ducks (Aix sponsa) inferred by the ‘isolation with migration’ coalescent method
.
Mol Ecol
 .
14
:
3407
3418
.
Ramos Pereira
MJ
Rocha
RG
Ferreira
E
Fonseca
C
.
2013
.
Structure of small mammal assemblages across flooded and unflooded gallery forests of the Amazonia-Cerrado ecotone
.
Biotropica
 .
45
:
489
496
.
Ramos-Onsins
SE
Rozas
J
.
2002
.
Statistical properties of new neutrality tests against population growth
.
Mol Biol Evol
 .
19
:
2092
2100
.
Ratter
JA
Bridgewater
SB
Ribeiro
JF
.
2003
.
Analysis of the floristic composition of the Brazilian Cerrado vegetation III: comparison of the woody vegetation of 376 areas
.
Edinb J Bot
 .
60
:
57
109
.
Rocha
RG
Ferreira
E
Fonseca
C
Justino
J
Leite
YL
Costa
LP
.
2014
.
Seasonal flooding regime and ecological traits influence genetic structure of two small rodents
.
Ecol Evol
 .
4
:
4598
4608
.
Rocha
RG
Ferreira
E
Oliveira
BMA
Martins
ICM
Leite
YLR
Costa
LP
Fonseca
C
.
2011
.
Small mammals of the mid-Araguaia river in Central Brazil, with the description of a new species of climbing rat
.
Zootaxa
 .
2789
:
1
34
.
Ronquist
F
Huelsenbeck
JP
.
2003
.
MrBayes 3: Bayesian phylogenetic inference under mixed models
.
Bioinformatics
 .
19
:
1572
1574
.
Smith
MF
Patton
JL
.
1993
.
The diversification of South American murid rodents: evidence from mitochondrial DNA sequence data for the akodontine tribe
.
Biol J Linn Soc
 .
50
:
149
177
.
Smith
TB
Saatchi
S
Graham
C
Salbbekoorn
H
Spider
G
.
2005
.
Putting process on the map: why ecotones are important for preserving biodiversity
. In:
Purvis
A
Gittleman
J
Brooks
T
, editors.
Phylogeny and conservation
 .
Cambridge
:
Cambridge University Press
. p.
166
197
.
Tajima
F
.
1989
.
Statistical method for testing the neutral mutation hypothesis by DNA polymorphism
.
Genetics
 .
123
:
585
595
.
Tamura
K
Stecher
G
Peterson
D
Filipski
A
Kumar
S
.
2013
.
MEGA6: molecular evolutionary genetics analysis version 6.0
.
Mol Biol Evol
 .
30
:
2725
2729
.
Valente
CR
Latrubesse
EM
.
2012
.
Fluvial archive of peculiar avulsive fluvial patterns in the largest Quaternary intracratonic basin of tropical South America: the Bananal Basin, central-Brazil
.
Palaeogeogr Palaeoclimatol Paleoecol
 .
356–357
:
62
74
.
Villela
FNJ
Nogueira
C
.
2011
.
Geologia e geomorfologia da estação ecológica Serra Geral do Tocantins
.
Biota Neotrop
 .
11
:
217
230
.
Voss
RS
Gutiérrez
EE
Solari
S
Rossi
RV
Jansa
SA
.
2014
.
Phylogenetic relationships of mouse opossums (Didelphidae, Marmosa) with revised subgeneric classification and notes on sympatric diversity
.
Am Mus Novit
 .
3817
:
1
27
.
Voss
RS
Jansa
SA
.
2009
.
Phylogenetic relationships and classification of didelphid marsupials, an extant radiation of new world metatherian mammals
.
B Am Mus Nat Hist
 .
322
:
1
177
.
Wallace
AR
.
1852
.
On the monkeys of the Amazon
.
Proc Zool Soc Lond
 .
20
:
107
110
.

Author notes

Corresponding Editor: Horacio Schneider