Amazonia as the Origin and Diversification Area of Didelphidae (Mammalia: Metatheria), and a Review of the Fossil Record of the Clade

Didelphidae is the largest New World radiation of marsupials, and is mostly represented by arboreal, small- to medium-sized taxa that inhabit tropical and/or subtropical forests. The group originated and remained isolated in South America for millions of years, until the formation of the Isthmus of Panama. In this study, we present the first reconstruction of the biogeographic history of Didelphidae including all major clades, based on parametric models and stratified analyses over time. We also compiled all the pre-Quaternary fossil records of the group, and contrasted these data to our biogeographic inferences, as well as to major environmental events that occurred in the South American Cenozoic. Our results indicate the relevance of Amazonia in the early diversification of Didelphidae, including the divergence of the major clades traditionally ranked as subfamilies and tribes. Cladogeneses in other areas started in the late Miocene, an interval of intense shifts, especially in the northern portion of Andes and Amazon Basin. Occupation of other areas continued through the Pliocene, but few were only colonized in Quaternary times. The comparison between the biogeographic inference and the fossil records highlights some further steps towards better understanding the spatiotemporal evolution of the clade. Finally, our results stress that the early history of didelphids is obscured by the lack of Paleogene fossils, which are still to be unearthed from low-latitude deposits of South America.


Introduction
The Metatheria clade corresponds to the extant marsupials plus all the extinct mammals more related to those than to the placentals (see Goin et al. 2016). Its oldest purported fossil (Sinodelphys szalayi, from the Early Cretaceous of Asia; Luo et al. 2003) was recently attributed to Eutheria, resulting in a 50-million-year ghost lineage (i.e., gap in the fossil record) for Metatheria (Bi et al. 2018). Metatherians formed a diverse group in Laurasia during the Late Cretaceous, where they are abundantly found, especially in fossil assemblages of North America. Currently, there is a single living group of Metatheria, the Marsupialia, which is taxonomically and morphologically diverse in Australia, South America, and Central America (Goin et al. 2016). The most likely hypothesis for the disappearance of metatherians in the northern hemisphere is related to the effects of the Cretaceous-Paleogene (K-Pg) extinction (Longrich et al. 2016).
The group probably reached South America following tectonic activities in the Caribbean plate that created a temporary bridge between South and North America during the Late Cretaceous (Campanian stage ;Pascual 2006). Case et al. (2005) considered that the most probable dispersal path was the Aves Ridge, a volcanic arc in the Caribbean Sea that began to dissipate in the Paleocene. Goin et al. (2012) named this migration episode as the First American Biotic Interchange (FABI), which allowed multiple dispersals of metatherians, among other vertebrates. Up to now, the oldest record of Metatheria in South America comes from the Paleogene of Argentine Patagonia (Cocatherium lefipanum; Goin et al. 2006).
Currently, there are over 130 species of marsupials in the Americas, and the group corresponds to around 10 % of the terrestrial mammal diversity in South America (Goin et al. 2016;Burgin et al. 2018). Most of them live in tropical and/or subtropical forests, are small-to medium-sized, nocturnal or crepuscular, and have arboreal habits; their offspring is very small and completely dependent on the mother, although only large-and medium-sized genera carry them in the marsupium, a pouch in the ventral region that gives the name to the clade (see Goin et al. 2016).
American marsupials are classified into three groups: Didelphimorphia, Paucituberculata, and Microbiotheria. The first is the most diverse, widely distributed in tropical and subtropical regions of South America, although one species (Didelphis virginiana) extends to southern Canada (Goin et al. 2016;Figueiredo and Grelle 2018). All living didelphimorphians, popularly known as opossums, are included in Didelphidae (sensu Voss and Jansa 2009). This clade probably originated in the Paleogene of South America and remained isolated in this continent for millions of years, until the formation of the Isthmus of Panama. The living Paucituberculata are grouped in only one family, Caenolestidae, with seven species attributed to Caenolestes, Rhyncholestes, and Lestoros, all distributed along the Andes (Goin et al. 2016). Conversely, extinct Paucituberculata were ecologically and taxonomically diverse during the Miocene of South America (Abello 2013). This is the same case for Microbiotheria, which contains a single living species, Dromiciops gliroides, restricted to the southern portion of the Andes (Martin 2010). This group has been subject of intense research in the last decades since Szalay (1982) proposed that Microbiotheria is part of the Australidelphia clade, indicating that it is more closely related to the Australian marsupials than to the South American ones, a hypothesis supported by several later studies based on molecular data (e.g., Kirsch et al. 1991;Mitchell et al. 2014).
Considering the relevance of marsupials in the South American fauna, the objectives of this work are: (i) to estimate the area of origin and diversification of Didelphidae based on parametric models for ancestral area estimates; (ii) to compile the pre-Quaternary fossil record of the group; and (iii) to contrast these data to each other and to major biogeographic events that occurred in South America during the Cenozoic.

Phylogenetic Tree
The time-calibrated tree topology adopted in our analyses is part of the phylogeny proposed by Mitchell et al. (2014), which presented a Bayesian inference of the relationships among extant marsupials, based on mitochondrial genomes and nuclear loci. It sampled 17 genera and 44 species of Didelphidae, which composed our ingroup; we selected the microbiotherian Dromiciops gliroides as the outgroup.

Biogeographic Areas and Distribution of Species
The biogeographic areas used herein were based on the ecoregions defined by Olson et al. (2001), the distribution pattern of didelphids, and the historical continuity of the landmasses from the Late Cretaceous to the Holocene, added to factors that may act as dispersal barriers for small mammals with limited swimming abilities. For example, even though the Didelphidae fauna of Central America and the lowlands northwestern to the Andes in South America are very similar, we consider those as two different units because they had no connection before the establishment of the Isthmus of Panama, at ca. 2.7 million years ago (Mya). Thus, we defined nine biogeographic areas along the current distribution of Didelphidae: (1) Central America and North America; (2) Northwestern South America, composed of trans-Andean humid-lowlands; (3) Andes, with montane formations above 1,500 m altitude; (4) Venezuela, composed of coastal forests and Llanos savannas; (5) Amazonia, composed of contiguous lowland rainforests alongside the Amazon, extending across the coastal drainages of the Guianas and northern Brazil; (6) Dry Diagonal, encompassing Chaco, Cerrado, Caatinga, and other seasonally dry tropical forests and savannas; (7) Atlantic Rainforest, lowland and montane forests along the Atlantic coast and the adjoining moist-subtropical forests of Paraná Basin; (8) Pampas, composed of the subtropical and temperate grasslands of southern Brazil and Uruguay; and (9) Steppes, temperate grasslands and shrublands, including Argentine Pampas and Patagonia (Fig. 1A). Fig. 1 Map of the biogeographic areas adopted in this study A, paleomaps of South America related to the episodes that defined the time bins of the analyses B, (generated using the softwares GPlates 2.0 and PaleoData Plotter; Scotese 2016), and time-calibrated tree of Didelphidae (Mitchell et al. 2014) with the ancestral areas reconstructions produced by our best-fit model C, Major subclades are indicated on the right, along with a photo of one of its taxa, boxes next to each terminal taxon represent its current area distribution, boxes on each node represent the most probable ancestral area, shades behind the tree indicate the stages of each epoch, and dashed lines separate the defined time bins (at 55, 10, and 2.7 Mya). All photographs are licensed under the Creative Commons Attribution-Share Alike 4.0 International license (see details in Supplementary File) ◂ 1 3 We established the presence of each species in the biogeographic areas by overlapping their present-day distributions to the areas map of Fig. 1A. Those distributions were based on the polygons available at the International Union for Conservation of Nature and Natural Resources website (IUCN 2021). We allowed a maximum range size of seven areas, based on the maximum joint distribution of extant taxa.

Ancestral Area Estimation and Biogeographic Events
To reconstruct the origin and diversification areas of Didelphidae, we used the package BioGeoBEARS (Matzke 2013) on the computing environment R (v3.5.1; R Core Team 2018). This package enables the evaluation of competing models for ancestral area reconstructions and stratified analyses over time. Two different basic models were employed using the same time-calibrated topology: DEC M0 (Dispersal-Extinction-Cladogenesis model; Ree 2005; Ree and Smith 2008) and DIVALIKE M0, an adaptation of the dispersal-vicariance analysis (DIVA; Ronquist 1997). Under both models ranges evolve through anagenetic stochastic changes (i.e., dispersals and extinctions, controlled by the parameters d and e, respectively) or three types of cladogenetic events, but DEC and DIVA differ regarding the implemented cladogenetic processes. DEC incorporates subset sympatry (one daughter species inheriting a subset of the ancestral area; e.g., ABCD -> A, ABCD) but disallows widespread vicariance (both daughter species inheriting joint ranges from a widespread ancestor; e.g., ABCD -> AB, CD), whereas in DIVALIKE widespread vicariance is available but subset sympatry is not (Ronquist and Sanmartín 2011;Matzke 2013). Narrow vicariance (daughter species inheriting non-overlapping subsets from a widespread ancestor; e.g., ABCD -> A, BCD) is allowed in both models. Additionally, the implementation of cladogenetic events based on biogeographic processes is the main difference between explicit biogeographic models (e.g., DEC, DIVA) and other types of analyses using methods not specifically designed to explore range evolution (e.g., ancestral state estimates using Maximum Likelihood).
We defined four time bins separated by three episodes that greatly changed the connectivity among the areas (Hoorn et al. 2010;Woodburne 2010;Goin et al. 2016;Scotese 2016;Fig. 1B). To determine the weight of transitions between the defined areas, we established multiplier matrices for each time bin, as follows: adjacent areas = 1.0; nonadjacent areas separated by only one other area = 0.5; adjacent areas separated by current (e.g., Andes and Paraná River) or past (e.g., marine transgression) geographic barriers = 0.5; areas separated by two or more areas = 0.
Time bin 1 (0-2.7 Mya): The relation between the areas is based on the current geography and distribution of didelphids. Its older boundary is marked by the establishment of the Isthmus of Panama, which resulted in the Great American Biotic Interchange (GABI), allowing the exchange of terrestrial faunas between North and South Americas.
Time bin 2 (2.7-10 Mya): the older boundary of this time slice is marked by the intense uplift of the northern Andes, causing profound environmental changes in the Amazon Basin. The relation between the areas is the same as in the first time bin, except for the unavailability of Central and North America and their connection to the northwestern South America area.
Time bin 3 (10-55 Mya): massive marine transgressions, caused by the sea level rise related to the Paleocene-Eocene Thermal Maximum (PETM), establish the beginning of this time bin. The Paranaense Sea enters the continent creating a barrier between the Steppes and both the Pampas and Atlantic Rainforest. The northern Andes is not uplifted yet, so there are no barriers between northwestern South America and both the Amazon and Venezuela areas.
Time bin 4 (55-95 Mya): before the PETM and the intense uplift of the Andean Cordillera, the Northwestern portion of South America is covered by a marine transgression and the connectivity between the Andes area and the adjacent areas (e.g., Steppes) is increased.
BioGeoBEARS allows two alternative ways to incorporate the weight of distance and connectivity: using symmetric matrices of absolute distances or distance multipliers between areas. The former is a novel feature of BioGeoBEARS and has been validated elsewhere (Van Dam and Matzke 2016). However, it requires the establishment of realistic distances between areas, which becomes increasingly difficult when we consider the geological changes of deep time. Distance multipliers, on the other hand -also available in other implementations of DEC (Matzke 2013) -influence the relative dispersal probability and, as such, are a proxy for the difficulty of a given organism changing from one area to another. The values chosen in the dispersal multiplier matrices are user defined. As there is no objective way to define the real weight of transition between two areas, they will inevitably incorporate some subjectivity. To evaluate the impact of our assumptions on these matrices, we employed two approaches: (1) inclusion of w as a free parameter (resulting in two additional models, DEC M1 and DIVALIKE M1) and (2) a set of sensitivity analyses. The w parameter is an exponent of the dispersal multiplier matrix (set to 1 as default) and, if inferred during the analysis, functions as a test on the weight of the matrix (Dupin et al. 2017;Marsola et al. 2019). For example, if w is inferred to be 0 then the matrix will have no effect on the ancestral area estimates. For the sensitivity analyses, we followed Poropat et al. (2016) in creating "harsh" and "relaxed" versions of our "standard" matrices, by respectively increasing and decreasing the weight of medium transitions (i.e., 0.5 values were transformed to 0 and 1, respectively). We also added two additional matrices to further explore their impact: "not-so-harsh" (0.5 transformed to 0.25) and "not-so-relaxed" (0.5 transformed to 0.75). The DEC + j model (Matzke 2013), a variant that emphasizes "jump dispersal" or founder-event speciation, was not tested here due to previous conceptual and statistical criticisms (Ree and Sanmartín 2018). For each of these diversification models, we calculated the Akaike Information Criterion (AIC) and the Likelihood Ratio Test, which allows the comparison and choice between competing models (Matzke 2014).

Origin and Early Diversification Areas of Didelphidae
DEC M0 (w = 1, standard matrix) and DIVALIKE M0 were the best fitted models, with DEC M0 yielding a slightly better AIC value (Table 1). Both models indicate the relevance of the Amazonia in the early diversification of Didelphidae. The most significant differences between these models were the inferred ancestral area of all didelphids and of the most recent common ancestor (MRCA) of Thylamys, Lestodelphys, Cryptonanus, and Gracilinanus. In both cases, DEC M0 estimated the Amazonia as the most likely ancestral areas and DIVALIKE M0 yielded larger joint areas, also including Amazonia (Supplementary Information). For such reasons we focus the rest of the analyses and discussions on the results of the DEC model (Fig. 1).
The estimates based on the sensitivity matrices generally agree with those of DEC M0. The exceptions were DEC M1 (+ w) and M2 ("harsh" matrix), which yielded extreme differences in relation to all other models. These results are coherent with the estimated value for w (= 1.221; see Supplementary Information), which upweighted the distance multiplier matrices, resulting in an effect similar to manually increasing the multipliers (i.e., similar to the expected effect of a "harsh" matrix). The probability of the best and alternative estimates are all similar and very low in most ancestral nodes, indicating that the uncertainty of the estimates was increased by more heavy distance multipliers. Noteworthy, the log-likelihood of both alternative models were significantly lower than that of the other scenarios (M1 = -245.75; M2 = -258.51), which is also indicated by the low-parsimonious inferences of the ancestral areas (Ree and Sanmartín 2018). Hence, the sensitivity analyses show our distance multiplier matrices have little impact on the results. DEC M0 (Fig. 1C) shows that the group originated and diversified within Amazonia for almost 30 million years. Occupation of other areas started only in the late Miocene. At ca. 10 Mya the MRCA of Didelphini expanded to a large joint area also encompassing northwestern South America, Andes, and Venezuela. Still during the late Miocene (ca. 7 Mya), the MRCA of the Thylamini Lestodelphys and Thylamys occupied the Steppes, whereas a subclade of Gracilinanus expanded its distribution to include the Andes. By the Miocene-Pliocene boundary (ca. 5 Mya), subclades of the Marmosini Marmosa expanded to the Andes, Venezuela, Atlantic Rainforest, and Dry Diagonal, whereas a subclade of the Caluromyinae Caluromys reached the Andes, Atlantic Rainforest, and Dry Diagonal areas, besides Amazonia in both cases. In the late Pliocene (ca. 3 Mya), the MRCA of the Thylamini Cryptonanus was present in Amazonia, Dry Diagonal, and Steppes, and a subclade of the Marmosini Monodelphis also occupied the Dry Diagonal. Cladogenetic events in Central and North America and in Pampas occurred only during the Quaternary.

The Fossil Record of Didelphidae
The present compilation aims to include all fossil records of Didelphidae available in the literature from the oldest, in the Miocene, until the end of the Pliocene; those are contrasted to the estimated ancestral areas in The Fossil Record of Didelphidae Versus Estimated Age and Ancestral Areas and in Fig. 2. The Quaternary records, although more numerous and widespread in South America, show a poorer temporal constraint; it means that most records are loosely dated as late Pleistocene or late Pleistocene-Holocene, lacking precise stratigraphic control during fossil collection and/or radioisotopic dates. It is worth mentioning that many of the fossil records correspond to dental remains and that very few were included in phylogenetic analyses so far; thus, the diagnostic features on which their taxonomic assignations are based have not been corroborated as synapomorphies yet. This review may be useful to future molecular divergence time analyses, as node calibrations are frequently based on fossils (e.g., Pavan  Dunn et al. 2013), with several isolated remains retrieved from the sediments of the Sarmiento Formation and tentatively attributed to the family (Goin et al. 2007;Beck and Taglioretti 2020).
The subsequent records come from the late Miocene of several localities. Antoine et al. (2016) mentioned a 'marmosid didelphimorphian' from an early late Miocene ("Mayoan") locality of the Pebas Formation, in the Peruvian Amazon.
From the Loma de Las Tapias Formation, San Juan Province, northwestern Argentina, the marmosin Thylatheridium sp. was retrieved as part of a Chasicoan assemblage, whereas an indeterminate didelphid was recovered in a Huayquerian assemblage (Contreras and Baraldo 2011).
Didelphid remains are also known from the late Miocene lower member of the Ituzaingó Formation, Mesopotamiense of Entre Ríos, Argentina. The age of this site is a matter of debate, but its maximum is ca. 9.5 Mya (Pérez 2013), whereas the mammalian fauna indicates a Huayquerian assignation (Cione et al. 2000). Goin et al. (2013) revised previous works and concluded that three didelphid taxa occur in this locality: Zygolestes paranensis (Thylamini), Philander entrerianus, and Chironectes sp. (both Didelphini).
With respect to late Miocene Amazonian records, Czaplewski (1996) described, but did not name, at least two didelphids based on isolated teeth from the Peruvian side of the Acre River ('Acre Conglomerate', Madre de Dios Formation; Huayquerian). For the same locality, Goin and de los Reyes (2011) described Lutreolina materdei, based on a single lower molar. From the Brazilian side of the Acre River, Cozzuol et al. (2006) described Didelphis solimonensis based on a partial mandible from the Solimões Formation (ca. 9-6.5 Mya; Huayquerian; Latrubesse et al. 2007).
Among the Thylamyni, Thylamys is found in the Montehermosan of Buenos Aires Province (Mones 1980;Deschamps and Tomassini 2016), and Lestodelphys in the Marplatan of the same province (Goin 1999;Martinelli et al. 2013); also, an indeterminate Thylamini is recorded in the late Pliocene (Marplatan) of the northwestern Jujuy Province (Ortiz et al. 2012).
The only possibly Pliocene record of Didelphidae outside Argentina comes from El Breal de Orocual, an asphalt deposit in northern Venezuela, where Didelphis sp. was retrieved from late Pliocene-early Pleistocene levels (?Marplatan; Rincón et al. 2009;Ruiz-Ramoni et al. 2018).
The following additional records should be taken into account if sparassocynids are considered to be part of Didelphidae, as indicated by the recent phylogenetic analysis of Beck and Taglioretti (2020), which places sparassocynids within or as sister to the extant Marmosini Monodelphis. Those carnivorous-adapted marsupials are classically considered non-didelphid didelphimorphians/didelphoids, usually with a family-rank (i.e., Sparassocynidae; see Beck and Taglioretti 2020). Indeterminate sparassocynid records include late Miocene remains from the Chasicoan of the Pampean Region (Goin 1995), and from the Huayquerian of Cerro Azul, La Pampa Province; also from this locality, Thylatheridium hudsoni may represent a sparassocynid, as mentioned above (Goin et al. 2000;Forasiepi et al. 2009). Hesperocynus is recorded in the late Miocene of the Cerro Azul Formation (La Pampa; Huayquerian), in the Andalhuala Formation (Catamarca Province; Huayquerian or Montehermosan-Chapadmalalan), as well as in the La Huertita Formation (Mendoza Province; Montehermosan-Chapadmalalan) (Forasiepi et al. 2009;Garrido et al. 2014;Chiesa et al. 2019). Apart from a species of questionable taxonomic status from the Montehermosan of Bolivia (Abello et al. 2015), Sparassocynus is recorded in the middle Miocene of Quebrada Honda, Bolivia (Laventan; Croft 2007), and in latest Miocene-Pliocene deposits of central and northwestern Argentina (Montehermosan, Chapadmalalan, and Marplatan; Forasiepi et al. 2009;Abello et al. 2015;Beck and Taglioretti 2020).

Origin and Early Diversification of Didelphidae: When and Where
The analysis of Mitchell et al. (2014) indicated that, after a long gap since Didelphimorphia diverged from other marsupials (in the Late Cretaceous), the first internal split of crown Didelphidae occurred in the late Eocene (ca. 38 Mya), and all major clades had diversified up to the middle Miocene. The results of Mitchell et al. (2014) roughly agree with others based on diverse molecular markers (e.g., Steiner et al. 2005;Meredith et al. 2011;Vilela et al. 2015), which range from the middle Eocene to the early Oligocene. On the other hand, Jansa et al. (2014), based on nuclear genes, inferred a younger basal split of Didelphidae, in the late Oligocene (ca. 26 Mya), similar to the age found by Beck and Taglioretti (2020;ca. 24 Mya). Previous phylogenetic analyses of crown Didelphidae coincide in the position of the major clades, except for the relationships between Caluromyinae and Glironiinae, which are either sister taxa or the former is the sister group of the remaining didelphids (e.g., Jansa 2009 contra Mitchell et al. 2014).
Most models evaluated here indicate the relevance of Amazonia in the early evolution of Didelphidae, hence we further discuss the most probable scenarios based on the DEC M0 model (Fig. 1C), as it provided the best statistical fit for our data. Out of the 43 cladogenetic events in the ingroup, 26 occurred in Amazonia, in addition to ten reconstructed to Amazonia along with one or more additional areas. Those cladogeneses include the origin of all major clades traditionally ranked as subfamilies (i.e., Caluromyinae, Glironiinae, Hyladelphinae, and Didelphinae). Within the most diverse of them (i.e., Didelphinae), the origin of the clades traditionally ranked as tribes (i.e., Marmosini, Thylamini, Metachirini, and Didelphini) ranged from the late Oligocene until the middle Miocene (i.e., ca. 26 to 10 Mya; Mitchell et al. 2014) and were all reconstructed as occurring in Amazonia as well, except for the origin of Didelphini, reconstructed to a joint area that, besides Amazonia, includes Northwestern South America, Andes, and Venezuela.
Our results indicate that Didelphidae began to diverge beyond Amazonia in the late Miocene, an interval of intense environmental changes in the northern portion of the Andes and in the Amazon Basin (Hoorn et al. 2010; see Environmental Events and the History of Didelphidae). During this interval, the didelphids occupied the northwestern South America, Andes, Venezuela, and Steppes areas. The earliest of these events, at ca. 10 Mya, corresponds to the origin of Didelphini in a compound area, as mentioned above.
Our results partly agree with the maximum likelihood reconstruction of habitats conducted by Jansa et al. (2014), which pointed out that the origin and early diversification of Didelphidae occurred in moist forests of South America, with subsequent dispersal to dry forests in the late early Miocene (ca. 16 Mya, by the lineage of the Marmosini Tlacuatzin), late middle Miocene (ca. 9 Mya, by the Thylamini Cryptonanus and Thylamys shortly after), and earliest Pliocene (ca. 5 Mya, by the Didelphini Lutreolina). Using Bayesian reconstruction of ancestral states, Mitchell et al. (2014) inferred that the origin and early diversification of Didelphidae occurred in wet-closed forests as well, with occupation of mesic-open environments becoming more likely from the early middle Miocene onwards (ca. 14 Mya). The above-mentioned studies, however, employed trait evolution reconstructions, which are not designed to address range evolution. In contrast, our analyses employ explicit biogeographic areas (i.e., nine areas; Fig. 1A), biogeographic models (i.e., DEC and DIVALIKE), and time-stratification (which allows modeling changes in the connectivity of the areas).
Recent estimations indicated that the MRCA of Didelphidae was quite small (22-33 g; Amador and Gianini 2016) and arboreal (see Goin et al. 2016: 60), which can represent an indirect support to an origin in forest environments. From a physiological perspective, extant didelphids generally do not tolerate cold-temperate environments, and small opossums (e.g., Thylamys and Lestodelphys) can inhabit higher latitudes/altitudes likely due to specific adaptations such as caudal fat storage and/or daily torpor (McNab 1986;Goin 1999). In this sense, it is interesting to note that didelphids reach their species richness climax at low latitudes of the Neotropics (Goin 1997(Goin , 1999Sánchez-Villagra et al. 2007), and the lineages Caluromyinae, Glironiinae, and Hyladelphinae are mostly restricted to tropical rainforests (Goin et al. 2016). Concordantly, a recent spatial analysis of species richness and phylogenetic diversity among Didelphidae indicated highest values in areas of tropical South American forests, whereas the lowest values were found at the northern and southern extremes of its distribution along the American continent (Figueiredo and Grelle 2018).
Regarding the short-tailed opossum Monodelphis, Vilela et al. (2015) inferred an early Oligocene (ca. 32 Mya) origin in the northern part of South America and Atlantic Forest. Pavan et al. (2016) also reconstructed a hybrid biome composed by Amazonia and Atlantic Forest as the area of origin of Monodelphis, but much later, in the late Miocene (ca. 7 Mya); the authors also concluded that these tropical lowland rainforests hosted the diversification of the genus until the Pleistocene, and were the primary sources for dispersal to other biomes. Despite the lower number of Monodelphis species sampled by Mitchel et al. (2014), our results indicate an Amazonian origin for the genus during the middle Miocene (ca. 12 Mya). Our analyses also show that the MRCA of Marmosa robinsoni and Ma. mexicana occupied the Venezuela area at the Miocene-Pliocene boundary (ca. 5 Mya).
With respect to Didelphis, Dias and Perini (2018) found that the basal-most divergence within the genus took place between the late Miocene and the early Pliocene (ca. 6 Mya), corresponding to the separation of D. virginiana from the Neotropical species. The authors claimed that this age is more consistent with the oldest fossil record of Didelphis (D. On the Thylamini, Jansa et al. (2014) inferred that its MRCA inhabited closed-canopy forests, and that colonization of open habitats on the east of the Andes occurred prior to the Lestodelphys-Thylamys split. Giarla and Jansa (2014) analyzed both genera, which primarily inhabit open biomes, and concluded that they originated in the lowlands eastern to the Andes during the Pliocene, later expanding into and across the Andes. They also found that rivers played an important role on speciation in lowlands, whereas differentiation among montane habitats may have affected speciation within the Andes, probably influenced by climatic events such as the Quaternary glacial cycles. Palma et al. (2014) also pointed out that rivers would be involved in the differentiation of Thylamys. The authors concluded that the Thylamys lineage diverged in the latest Oligocene (ca. 24 Mya) and had the first internal split during the early Miocene (ca. 19 Mya) in the Chacoan ecoregion (sensu Morrone 2006, which roughly encompasses the areas defined herein as Dry Diagonal, Pampas, and northeastern part of the Steppes). The analyses of Vilela et al. (2015) also indicated that most divergences within Thylamini occurred in the Miocene, with the exception of the Pliocene diversification of Cryptonanus. Our results agree with the inference of Jansa et al. (2014) that Thylamyni originated in closed forests, and estimate Amazonia as their ancestral area; also, the colonization of open biomes before Lestodelphys-Thylamys split is corroborated. Our analyses indicate that the MRCA of both genera was distributed in the Steppes in the late Miocene (ca. 7 Mya). On the diversification of Thylamys, despite the age estimations of Palma et al. (2014) being considerably older (with large uncertainty) than those of Mitchell et al. (2014, which indicate the first internal split of the genus in the late Pliocene), their biogeographic reconstruction of a Chacoan origin partly agrees with ours, which points the origin of Thylamys in the Andes, Dry Diagonal, and Steppes area. Our results also indicate that, in the late Miocene (ca. 7 Mya), a subclade of Gracilinanus that includes G. agilis, G. microtarsus, and G. aceramarcae expanded its estimated ancestral area to include Amazonia and Andes.

The Fossil Record of Didelphidae Versus Estimated Age and Ancestral Areas
The comparison of each fossil record of Didelphidae detailed in The Fossil Record of Didelphidae Versus Estimated Age and Ancestral Areas with the estimated age of its clade (Mitchell et al. 2014) and with the ancestral areas reconstructed by our best-fit model (DEC M0 model) is given below, organized by the age of the fossils, and illustrated on Fig. 2. Sparassocynids were recently included in Marmosini (Beck and Taglioretti 2020), but this systematic proposition has not been scrutinized and adopted by the specialists yet. Therefore, we discuss these taxa as tentative marmosines, differentiating their records (i.e., dashed bones in Fig. 2) and discussing them at the end of this section. In our comparison, four scenarios were identified: (1) the fossil record agrees with both the age and the ancestral range estimated for a clade; (2) the fossil record is older than the age estimated for a clade; (3) the fossil record indicates an earlier occupation of an area than our biogeographic model; or (4) the fossil record documents the presence of a clade in an area where it is currently absent. It is clear that incongruences between the fossil record and the inferred ages/ancestral ranges can be due to multiple causes, including mistaken taxonomic assignation and inaccurate temporal/spatial estimations. Anyway, our results stress the importance of a timetree that encompasses the fossil record of Didelphidae in order to better understand the spatiotemporal evolution of the clade.
Our biogeographic reconstruction shows that Didelphidae was restricted to low latitude humid forests of Amazonia since its origin, in the late Eocene, up to the late Miocene, when different lineages started to diverge in other areas. The oldest putative record of Didelphidae comes from the early Miocene (Colhuehuapian, ca. 21.1-20.1 Mya) of Argentine Patagonia (Goin et al. 2007), which is encompassed in the Steppes of our analysis. However, our estimates indicate that this area was colonized only in the late Miocene (ca. 7 Mya), thus representing an example of scenario 3. The wide temporal and spatial gap suggests that the record from the Colhuehuapian sediments of Gran Barranca represent either a Didelphoidea, as suggested by the original authors (Goin et al. 2007), or a lineage of crown Didelphidae that left no extant representatives. The absence of didelphids in subsequent early and middle Miocene Patagonian assemblages (including the mammal-rich Santacrucian) seems to favor the first hypothesis.
Contrasting the oldest definitive records, from the middle Miocene (15.6-11.8 Mya;Marshall 1976;Goin 1997;Antoine et al. 2013), and the estimated age for the first internal split of the crown clade (38 Mya, Mitchell et al. 2014), the ghost lineage of Didelphidae is ca. 24 million years. If we consider alternative age estimates for the same event, which range from the middle Eocene to the late Oligocene Vilela et al. 2015), this ghost lineage spans from 30 to 12 million years. In any case, the early evolutionary history of didelphids is obscured, as Paleogene fossils are unknown. This probably relates to the paucity of Eocene and Oligocene localities at low latitudes of South America, which is the estimated area of origin and early diversification of Didelphidae in our analyses, as well as in others. Previous authors also recognized that the relatively poor Paleogene fossil record of tropical South America hinders understanding the evolutionary and biogeographic history of the clade (Cozzuol et al. 2006;Sánchez-Villagra et al. 2007;Antonelli et al. 2018;Beck and Taglioretti 2020).
These oldest definitive fossils were retrieved from La Venta, Colombia (Marshall 1976;Goin 1997) and Madre de Dios, Peru (Antoine et al. 2013), which correspond respectively to the areas of Andes and Amazonia. Both sites record Marmosa laventica; our analyses indicate that the MRCA of Marmosa occurred in Amazonia during the middle Miocene (ca. 15 Mya). Therefore, the earliest fossil record of Didelphidae in the Peruvian Amazonia agrees with our biogeographic reconstructions, exemplifying scenario 1. On the other hand, the record in the Andes predates the occupation of this area by the Marmosini, estimated here to the Miocene-Pliocene boundary (ca. 5 Mya), representing scenario 3; however, it corresponds to the time lapse in which the area could have been colonized by the ancestors of Marmosa (dashed segment in Fig. 2).
The record of Thylamys (i.e., Thylamys minutus and Thylamys colombianus) in La Venta exemplifies scenario 3, as the occupation of the Andes by the genus was estimated to be much later, in the late Pliocene (ca. 3 Mya). Amazonia was estimated as the ancestral area for Thylamini (in the Oligocene-Miocene boundary) and for the clade that includes Thylamys, Lestodelphys, Cryptonanus, and Gracilinanus (in the middle Miocene). Later, in the late Miocene (ca. 7 Mya), the Thylamys + Lestodelphys lineage became restricted to the Steppes. Beck and Taglioretti (2020) also pointed the temporal incongruence between this middle Miocene record of Thylamys and its divergence time (estimated as latest Miocene or Pliocene by their study), and suggested that it may not belong to the genus.
With respect to late Miocene records, the presence of a 'marmosid' didelphid in the Peruvian Amazon (Antoine et al. 2016) agrees with the estimated age and ancestral area of Marmosini, as by this time this lineage was restricted to the Amazonia (scenario 1). The disputed mandibular ramus from the Chasicoan of central Argentina (Reig 1957) suggests scenario 3, considering that the Steppes area was first colonized later, at ca. 6 Mya (i.e., Huayquerian-Montehermosan boundary, by the MRCA of Thylamys + Lestodelphys); however, the record falls within the time lapse when this area could have been colonized, as it postdates the preceding node (i.e., the split between the lineages Thylamys + Lestodelphys and Gracilinanus + Cryptonanus in the Amazonia at ca. 15 Mya). The Chasicoan record of the extinct Marmosini Thylatheridium in the Andes exemplifies scenario 3, as the occupation of the area was inferred to the Miocene-Pliocene boundary, by a subclade of Marmosa, and the phylogenetic placement of Thylatheridium is still to be tested.
Among the didelphids recovered from the Mesopotamiense of Argentina, the extinct Thylamini Zygolestes is classified as closely related to Gracilinanus (Goin et al. 2013), although this was not tested by a phylogenetic analysis. By this time (late Miocene, Huayquerian; Cione et al. 2000), the lineage of Gracilinanus was diversifying in Amazonia, whereas the fossil site is located in the Steppes area, posing a spatial incongruence (scenario 3). On the record of Chironectes and an extinct species of Philander, the origin of crown Didelphini (the clade encompassing both) is estimated to a joint area including northwestern South America, Andes, Venezuela, and Amazonia, slightly before that (at ca. 10 Mya). However, the estimated divergence of Philander is younger, in the Pliocene (ca. 4 Mya), posing a temporal incongruence for this record (scenario 2); Beck and Taglioretti (2020) also noticed that this fossil record of Philander predates the divergence time they estimated for the genus. In terms of space, extant species of the genus do not reach the Steppes area, same case of extant Chironectes (IUCN 2021), configuring a scenario 4. With regard to the Huayquerian didelphids of Cerro Azul, in central Argentina (Goin et al. 2000), the record of the Thylamini Thylamys is also older than the estimated origin of the genus (scenario 2). The record of the extinct Didelphini Hyperdidelphis is potentially in agreement with the estimated age of the tribe (at ca. 10 Mya), but would represent a much earlier invasion of the Steppes by the clade (scenario 3); in the case of Lutreolina, this Huayquerian fossil agrees with the estimated age of the genus (ca. 7 Mya), but not with the ancestral estimated area (Venezuela). On the other hand, the late Miocene Huayquerian record of an extinct species of Lutreolina in the Peruvian Amazonia (Goin and de los Reyes 2011) matches the spatiotemporal inference (scenario 1). The Huayquerian record of an extinct species of Didelphis in the Brazilian Amazonia (Cozzuol et al. 2006) slightly predates the estimated divergence time of D. virginiana (ca. 6 Mya); Beck and Taglioretti (2020) noticed the same regarding this last fossil record and the divergence time of the genus estimated by them.
On the Pliocene records, the presence of Marmosa and the extinct genus Thylatheridium in the Chapadmalan of Buenos Aires, Argentina, incurs in scenario 4, as sampled Marmosini do not occupy the Steppes area (but see below). Among the Thylamyni, the record of Thylamys and Lestodelphys in the Montehermosan and Marplatan, respectively, of Buenos Aires Province is in agreement with our estimates, which indicate that the MRCA of both genera occupied the Steppes earlier (at ca. 7 Mya; scenario 1); the same applies to the indeterminate Thylamini from the Marplatan of the Dry Diagional in Jujuy Province. Among the Didelphini, the presence of Lutreolina in the Montehermosan and Chapadmalalan of central Argentina agrees in the temporal aspect, as the genus diverged earlier than that, and the extant species occupy the Steppes area (scenario 1). The presence of Didelphis in the Chapadmalalan and Marplatan of central Argentina agrees with the divergence age of sampled South American species of genus (in the late Pliocene, ca. 3 Mya), and slightly predates the estimated occupation of the Steppes (scenario 3). The fossils of Didelphis retrieved from the late Pliocene-early Pleistocene of Venezuela are congruent with the estimated divergence of the South American species of the genus, as well as with the estimated occupation of the area (scenario 1). The evaluation of the records of the extinct genera Thylophorops and Hyperdidelphys in the Montehermosan and Chapadmalalan of central Argentina depend upon their phylogenetic placement, but also predates the estimated occupation of the area by the Didelphini.
Sparassocynids were recently placed within Marmosini by the phylogenetic analyses of Beck and Taglioretti (2020), as part of Monodelphis or as a sister group of this genus. Their fossil records are restricted to two areas of our biogeographic scheme. In the Andes, they correspond to Sparassocynus in the middle Miocene (Laventan; Croft 2007) and Miocene-Pliocene boundary (Abello et al. 2015), and Hesperocynus in the Huayquerian or Montehermosan-Chapadmalalan (Forasiepi et al. 2009;Garrido et al. 2014). Our analyses estimate the arrival of Marmosini to the Andes by the Miocene-Pliocene boundary, by a subclade of Marmosa. The sampled species of Monodelphis are absent in the Andes, but few other extant species of the genus occur peripherally in this area (IUCN 2021). In the Steppes, the records correspond to two indeterminate sparassocynids from the Chasicoan (Goin 1995) and Huayquerian (Goin et al. 2000), Hesperocynus from the Huayquerian (Forasiepi et al. 2009) and Montehermosan-Chapadmalalan (Garrido et al. 2014;Chiesa et al. 2019), and Sparassocynus from the Montehermosan-Marplatan (Forasiepi et al. 2009;Abello et al. 2015). Those records, as well as that of Marmosa and Thylatheridium in the Montehermosan-Marplatan (Marshall 1976;Goin 1999;Deschamps and Tomassini 2016), support the presence of Marmosini in the Steppes during the late Miocene and Pliocene (see below).
To sum up, some fossil records predate the estimated occupation of areas (scenario 3; records 3, 4, 7, 9-12, and 21-24 in Fig. 2). However, most of these (records 3, 7, 9, 11, and 21-24) are within the time lapse between the estimated ancestral occupation of the area and its previous node (dashed rectangles in Fig. 2), corresponding to the uncertainty (i.e., ancestral branch) within which the dispersal should have happened according our hypothesis. On the other hand, the remaining records (i.e., 4, 10, and 12; Fig. 2) represent "harder" incongruences, that might be related to taxonomic misinterpretations or to underestimation of the divergence dates due to the incomplete specieslevel sampling and/or the lack of fossils in the adopted timetree (Mitchell et al. 2014). Two other records are placed in areas where the sampled extant taxa do not occur (scenario 4; records 19 and 20), suggesting the past presence of Marmosini in the Steppes. Nevertheless, Monodelphis dimidiata (not included in the analysis of Mitchell et al. 2014) occur in the Steppes of Argentina (IUCN 2021), indicating the potential ability of the clade to thrive in this colder/drier environment; this is also indicated by the presence of sparassocynids in the Steppes and Andes areas (records 26-33 on Fig. 2). Regarding the records of indeterminate Didelphidae (records 1, 5, 6, 8, and 15 on Fig. 2), with the exception of the first one (discussed in the beginning of this section), all are younger than the estimated arrival of at least one subclade in their areas. It is also noticeable that Caluromyinae is estimated to occur in the Amazonia since the late Eocene, and also in the Andes, Atlantic Rainforest, and Dry Diagonal at least since the Miocene-Pliocene boundary; however, Caluromyinae has no known fossils, suggesting that there are factors reducing the fossilization potential of the group. In the same sense, the occupation of the Atlantic Rainforest by didelphids is estimated to the Miocene-Pliocene boundary, but the area lacks fossil sites that document this time lapse.

Environmental Events and the History of Didelphidae
In order to understand how the remarkable Neotropical biodiversity was assembled, many authors have analyzed the relative importance of Quaternary glacial/interglacial cycles versus older Paleogene and/or Neogene geomorphological events (e.g., closure of Panama Isthmus, Andean uplift, marine incursions; Rull 2011). In this sense, after 1 3 the Pleistocene refuge hypothesis was intensely contrasted to geological, paleontological, and molecular data, recent reviews converge that the diversification in the Neotropics was a complex process, shaped both by Quaternary climatic cycles and older geological events (Antonelli et al. , 2018Rull 2011;Werneck et al. 2012).
Previous research already pointed out that didelphids were affected by environmental changes, such as the uplift of mountain chains and repeated expansions-contractions of tropical forests and open formations (Costa 2003;Giarla and Jansa 2014;Figueiredo and Grelle 2018). In this section, we summarize the major Cenozoic environmental events that might have impacted the spatiotemporal diversification of Didelphidae.
Recently, Antonelli et al. (2018) showed that Amazonia was the region that most contributed to the total biodiversity of tropical America during the Cenozoic, as it not only generated enormous diversity in situ, but also provided lineages to all other Neotropical regions. Also, Rangel et al. (2018) showed that, under the influence of the Andes, western Amazonia was a key region to South American biodiversity. A review of the fossil record of western Amazonia also indicated that it was a hotspot of mammalian biodiversity throughout the Cenozoic (Antoine et al. 2017). By the estimated time of Didelphidae origin (late Eocene, ca. 38 Mya), a rainforest already covered the Amazonia for a long time, at least since the Paleocene (Wing et al. 2009;Woodburne et al. 2014). Also, it was part of a larger "pan-Amazonian" region, which included the present Amazon, Orinoco, and Magdalena basins, with some parts alternating between fluvial and marginal marine conditions (Hoorn et al. 2010). In the latest Oligocene, this forested lowland began to experience the Andean uplift (Hoorn et al. 2010). Climatic fluctuations also took place since the Paleogene-Neogene boundary, when the warm and wet forests started to give place to savannas (Goin et al. 2016). Previous authors already questioned if those climatic and vegetation shifts triggered the occupation of new areas by the Didelphidae (Steiner et al. 2005;Goin et al. 2007).
Our analyses emphasize the importance of the Miocene in the diversification and biogeographic dynamics of Didelphidae. For instance, almost half of the cladogenetic events sampled in the timetree occurred during this epoch, including the origin of several extant genera. In this sense, Rull (2011) pointed out that most Neotropical genera initiated their diversification in the Neogene and attained the current diversity in the Quaternary. Our biogeographic analyses indicate that cladogeneses in new areas started in the late Miocene, in agreement with the results of Antonelli et al. (2018), who found an increased exchange of animal and plant taxa between the Amazonia and other regions since this period. Under the influence of the Middle Miocene Climatic Optimum (MMCO; ca. 14-17 Mya) and the wetlands of the "Pebas System" in western Amazonia (Hoorn et al. 2010), the oldest definitive fossils of Didelphidae are inferred as occupying a mixture of tropical rainforest and open habitats under a monsoonal-like tropical climate (Antoine et al. 2013).
Intense environmental changes took place after the MMCO, such as the global sea level drop and cooling process, most intense peaks of Andean building (at ca. 12 Mya), establishment of the "Acre System" in western Amazonia, and the Amazon drainage to the Atlantic Ocean (at ca. 10 Mya; Hoorn et al. 2010). Concomitantly, the spread of open habitats in South America culminated in a clearer latitudinal zonation and the establishment of the modern version of shrublands and grasslands (see Goin et al. 2016). At the late Miocene, fragmented distributions induced by climate cooling in tropical, subtropical, and mesic habitats likely triggered the final burst of diversification in some taxa, e.g., modern birds (Claramunt and Cracraft 2015). Also, the fossil record shows a moment of high diversification of mammals in tropical lowlands by this time (ca. 11.5 Mya; Hoorn et al. 2010). In the late Miocene, the continuous uplift of the northern-most part of the Andes and the expansion of alpine environments (Hoorn et al. 2010) may have caused vicariant events. Also during this time, the Dry Diagonal experienced the final uplift of the Central Brazilian Plateau (ca. 7-5 Mya), followed during the Pliocene by the expansion of grass-dominated habitats and the increase of fire frequency (at ca. 5-3 Mya; Werneck 2011).
In higher latitudes, the late Pliocene shows a progressive cooling, also evidenced by floristic and faunistic changes, such as the appearance of the first steppe-adapted mammals (Tonni et al. 1992). According to Goin (1999), the fossil record of marsupials in central Argentina agrees with a faunal turnover postulated for South American terrestrial mammals close to the Pliocene-Pleistocene boundary (ca. 2.6-2.4 Mya), likely related to cooler and dryer conditions; for instance, it shows the last record of specialized, largesized carnivorous taxa (e.g., Thylophorops, Hyperdidelphys, sparassocynids).
Finally, the Quaternary glacial/interglacial cycles generated retractions/expansions of humid forests over savannas, which possibly correlates with cladogenetic events and 1 3 occupation of new areas by the Didelphidae (Goin et al. 2016;Dias and Perini 2018), as suggested for other taxonomic groups (see Rangel et al. 2018). Also, the influence of North American immigrating fauna is more prominent since the Plio-Pleistocene boundary (Tonni et al. 1992). Although most cladogenetic events in the adopted timetree are pre-Pleistocene, its modest taxonomic sampling at species-level precludes addressing the impact of the Quaternary cycles in the spatiotemporal evolution of Didelphidae (for this purpose, see e.g., Vilela et al. 2015 andPavan et al. 2016 for Monodelphis; Dias and Perini 2018 for Didelphis; Giarla andPalma et al. 2014 for Thylamys).

Final Remarks
Besides being the world's most diverse rainforest, Amazonia is proving to have been an important source of lineages to other Neotropical regions (Antonelli et al. 2018;Rangel et al. 2018). Our results show that this was also the case for Didelphidae, which originated and diversified in Amazonia for at least the first half of its history, before reaching other South American areas. Further studies aiming to improve the understanding of the spatiotemporal dynamics of Didelphidae should test if current species diversity and endemism are related to areas of historic stability or to heterogeneous areas, and analyze the impact of Quaternary climatic cycles. In a larger temporal scale, it would be valuable to conduct biogeographic analyses that focus on Amazonia, splitting it into smaller areas and testing the role of Cenozoic tectonic and paleogeographic events (e.g., Andes uplift, establishment of Pebas and Acre systems, changes in the drainage patterns). Another important advance would be the construction of a more inclusive timetree. This effort should improve the sampling of extant species, especially in light of the increased number of recognized taxa and increased availability of molecular data in recent years. Also, the taxonomic assignment of fossil taxa should be tested within a phylogenetic frame, ultimately building a timetree that encompasses extinct taxa as well. The complete compilation of the Didelphidae fossil record presented here can contribute to this, in addition to providing calibration dates to future molecular divergence time analyses. Finally, our analysis also highlights the lack of fossils that document the initial history of Didelphidae, which are still to be unearthed from low-latitude Paleogene and early Neogene deposits of South America.
Funding Open Access funding enabled and organized by Projekt DEAL. No funding was received for conducting this research.
Data Availability All datasets analyzed during this research are available in the Supplementary Information file and on the GitHub repository, at https:// github. com/ gsfer reira bio/ didel phidae_ biogeo.

Declarations
Competing Interests The authors have no conflicts of interests related to this research.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.