Species co-occur with different sets of other species across their geographical distribution, which can be either closely or distantly related. Such co-occurrence patterns and their phylogenetic structure within individual species ranges represent what we call the species phylogenetic fields (PFs). These PFs allow investigation of the role of historical processes—speciation, extinction and dispersal—in shaping species co-occurrence patterns, in both extinct and extant species. Here, we investigate PFs of large mammalian species during the last 3 Myr, and how these correlate with trends in diversification rates. Using the fossil record, we evaluate species' distributional and co-occurrence patterns along with their phylogenetic structure. We apply a novel Bayesian framework on fossil occurrences to estimate diversification rates through time. Our findings highlight the effect of evolutionary processes and past climatic changes on species' distributions and co-occurrences. From the Late Pliocene to the Recent, mammal species seem to have responded in an individualistic manner to climate changes and diversification dynamics, co-occurring with different sets of species from different lineages across their geographical ranges. These findings stress the difficulty of forecasting potential effects of future climate changes on biodiversity.
Large-scale patterns of species diversity result from the interplay between ecological and evolutionary processes acting at different spatial and temporal scales [1,2]. Ultimately, historical factors are responsible for the richness and composition of regional biotas from which local communities form [3,4]. Consequently, there is an increasing interest on evolutionary approaches to understand biodiversity [4–6]. The evolutionary component of biodiversity patterns have been traditionally assessed by reconstructing the phylogenetic relationships among extant species or by using the fossil record . Relying on one type of historical information only (i.e. molecular phylogenies or fossil occurrences data) biases our understanding of evolutionary influences on biodiversity. Thus, recent claims favour a formal integration of the neontological and palaeontological perspectives to fully understand the patterns and processes in biodiversity over time and space [7–11]. In fact, several studies of current biodiversity now routinely include fossils, or apply phylogenetic comparative methods directly to fossil data (see  and references therein). Most of these studies, however, mainly focus on the temporal dimension of biodiversity patterns, neglecting the spatial dimension that sets the physical stage for evolutionary processes to take place .
The analysis of the spatial dimension in palaeontological studies is considered as the ‘final frontier’ [12,13], given that most of these studies involve local communities or global compendia of fossil data, without focusing on intermediate scales and often ignoring space altogether. Spatially explicit analyses of biodiversity patterns are one of the main focuses of macroecology, which is increasingly including species' evolutionary dynamics, either in the form of phylogenies (‘phylogenetic macroecology’) or fossil data (‘palaeomacroecology’; sensu ). A palaeomacroecological approach is thus a way towards that ‘final frontier’. One way of putting together palaeontological information and methods into a temporally and spatially explicit framework is to apply the burgeoning field of community phylogenetics  to study the assembly of fossil communities (e.g. ). This approach requires considering the phylogenetic structure of local communities relative to the species pool where they belong . However, community phylogenetics works with individual sites or communities (i.e. the community-level) that may be inappropriate to reveal species-level patterns unless all sites/communities where species occur (i.e. their geographical ranges) are considered .
Species' geographical ranges vary widely in size, shape and location, and their overlap ultimately determines the geographical variation in species richness at the macroecological scale . Thus, combining the spatial and temporal perspectives by considering the phylogenetic structure of species co-occurrences helps in evaluating the importance of the evolutionary process in driving geographical patterns of biodiversity [13,16]. Following this reasoning and relying on the community phylogenetics approach , Villalobos et al.  recently introduced the concept of ‘phylogenetic fields’ (PFs) to study patterns of species co-occurrence under a phylogenetic perspective. Instead of focusing on particular sites or regions as in traditional community phylogenetics, the PF concept focuses on the complete geographical ranges of species by describing the phylogenetic structure of species co-occurrences within individual ranges [16,18]. Thus, the PF is a species' property depicting the community, in terms of composition and history, with which a given species co-occurs . Consequently, PFs link the spatial and temporal dimensions to study broad-scale biodiversity patterns, allowing the investigation of evolutionary processes such as speciation, extinction and dispersal under appropriate macroecological and macroevolutionary scales [16,18]. However, species' ranges are not static. They usually show considerable temporal dynamics involving shifts, contractions and expansions [19–23]. Therefore, by including spatially explicit fossil data, it becomes necessary to study the temporal dynamics of PFs and to determine the role of historical and evolutionary processes in shaping geographical patterns of species co-occurrences.
Range shifts have been intensively studied owing to their relationship with historical and recent climate changes affecting biodiversity . Species' responses to such climate changes are fundamental to understanding the evolution of species and their organization into communities over time [22,25]. The prevailing idea is that species respond in an individualistic manner to changes in climate [26,27], generating great variability in community composition over time. For instance, in mammals, which are considered a model system for palaeontological studies owing to their relatively high preservation potential and commonness , extensive faunal turnover has been documented in different regions, especially over periods of intense climate changes such as the Quaternary . Despite high faunal turnover and individualistic responses of mammal species, including phylogenetic information into palaeocommunity analyses has revealed some interesting patterns. For example, Raia  showed a tendency of Eurasian Ice-Age mammal communities to become phylogenetically clustered towards the Recent, with particular clades being significantly overrepresented, suggesting that clade-by-clade replacement, rather than simple species turnover, took place . This and other studies focused either on particular sites and communities or on species within specific regions (e.g. North America; [25,26]), potentially missing part of species' ranges. Therefore, it remains to be determined if species' ranges and the phylogenetic structure of their co-occurrence patterns have also responded in an individualistic manner or if there were any congruent patterns among species.
Here, we apply a palaeomacroecological approach to study the temporal dynamics of co-occurrence patterns among mammalian species from the Late Pliocene to the Recent. We estimated PFs of extinct, large mammal species (body mass more than 5 kg) since 3 Ma to the Recent, and compare their temporal dynamics against diversification rates through the same time period, as estimated directly from fossil occurrence data. Based on previous findings of mammalian species and communities over time, we can expect two alternative scenarios to be true: (i) either idiosyncratic patterns of species' PFs, indicating individualistic responses, or (ii) a trend for clustered PFs towards the Recent species ranges, especially around times of dramatic climate changes, which would agree with clade-by-clade replacement, suggesting that habitat filtering is more likely to have occurred.
2. Material and methods
(a) Species occurrences
We assembled a dataset including 12 867 fossil occurrences for 313 Eutherian mammal species, inclusive of extinct taxa and some living species geographically spanning over the Palearctic, Afrotropic, Indomalaya, Nearctic and Neotropic regions. This dataset ranges in time from approximately 3 Ma to the Recent and includes minimum and maximum ages for each occurrence point per species. Following common palaeontological practice, we focused on mammal species with higher preservation potential, including species with body mass higher than 5 kg and excluding those in the orders Rodentia and Chiroptera as well as other clades restricted to forested habitats where fossilization is rare (e.g. Primates). Thus, as in most palaeontological studies of mammals from the Pliocene to the Recent, we concentrated on species from the orders Carnivora, Cetartiodactyla, Perissodactyla and Proboscidea. Such a criterion guarantees the consideration of species with comparable preservation potential, meaning that the fossilization bias is evenly distributed among taxa .
We collected information about fossil localities and their faunal lists from the following databases: http://paleobiodb.org (mainly for American localities), NOW database (http://www.helsinki.fi/science/now/, mainly for Eurasian and African localities), Miomap (http://www.ucmp.berkeley.edu/miomap/, for North American localities) and Jousse & Escarguel  for Holocene African localities. Although many absolute age estimates are available, when dealing with age-missing localities we used Spectral Ordering age estimation , an analytic method based on the similarity of localities' faunal lists [23,31].
Despite including both extinct (193) and living species (120) in our data, we only focused on extinct species for our co-occurrence and phylogenetic structure analyses (see below) to avoid the potential confounding factor of living species having their distributional patterns affected by anthropogenic activities rather than natural phenomena. Furthermore, species have predictable occupancy and geographical range temporal courses [23,31,32] and the inclusion of both extinct and extant species in the same analyses could yield results biased by the incompleteness of living species' evolutionary trends.
(b) Phylogenetic data
We obtained a phylogeny with 126 species included in the orders Carnivora, Cetartiodactyla, Perissodactyla and Proboscidea, for which we had at least three occurrences in the database (since it is not possible to estimate range size with less than three occurrences). The phylogeny was built upon the tree provided by Raia et al. , whose topology and (calibrated) ages of higher level clades were in turn derived from a species-level mammal supertree . This supertree combines different phylogenetic hypotheses (molecular and morphological) into a single tree, calibrated with both molecular and fossil dates (see  for details). Ages of inner nodes not taken from Bininda-Emonds et al.  were retrieved directly from the fossil record . Individual sub-phylogenies of fossil mammals were subsequently added according to specialists' opinions (see  for details). In our phylogeny, tips end at the last known appearance date of the species in the fossil record. Some intrafamilial relationships are not known and were therefore left as soft polytomies. This is known to have minor impact on variance–covariance matrix estimation [34–36], which is important for phylogenetic comparative methods . Species first appearance records in the NOW and PaleoDB databases were used to calculate branch lengths [34–36].
(c) Species diversity and distribution across time
To evaluate the temporal dynamics of species co-occurrence and its phylogenetic structure patterns, we first partitioned the fossil record in 0.5 Myr long consecutive time bins resulting in six time bins from 3 Ma to the Recent. Then, we selected all species present in each time bin and generated a phylogeny for each time bin by pruning the initial tree of 126 species (see above) to include only the species within each bin. Similarly, for each time bin and using the geographical coordinates of species occurrences, we generated a presence–absence matrix (PAM) representing species diversity (richness) and distribution (geographical ranges) across localities. To do so, we overlaid a regular grid of 1° resolution onto the occurrence points, which is a common procedure and resolution applied in macroecological studies . With this approach, we minimized the number of occurrences per species and simultaneously maximized the potential co-occurrence among species within grid cells. We considered point occurrences directly rather than creating minimum convex polygons (MCPs) to describe species' geographical distributions. Although MCPs are now routinely used in palaebiological studies, they inevitably overestimate species' geographical ranges, especially when their distribution is discontinuous, as frequently observed when analysing fossil occurrences [13,23,38].
From the PAM of each time bin, we described species' distribution and co-occurrence patterns. First, we estimated the species' range sizes per time bin. Then, we described species' co-occurrence patterns in two complementary ways: (i) by calculating the species richness comprised within individual species' ranges (hereafter, within-range richness) via recording the total number and identity of species occurring across all cells occupied by each species (species' diversity fields; ) and (ii) by estimating the C-score  between all species pairs, which is a commonly used index to quantify species association . We depicted the temporal patterns of range size, within-range richness and C-scores across time bins with box plots, allowing us to present variability of species within each time bin. In addition, to statistically compare such temporal patterns, we conducted Kruskal–Wallis rank sum tests followed by non-parametric multiple comparisons among time bins based on global ranks and Tukey's contrasts . In addition, we estimated the average proportional richness from the PAM of each time bin, which equals the reciprocal of Whittaker's beta diversity  thus providing a metric of species turnover.
Although we based our analyses on species with at least three occurrences in our database, we repeated the aforementioned analyses considering only those species with seven or more occurrences. This higher cut-off (seven occurrences) empirically represents a better compromise between small and large range sizes and sample size within time bins (i.e. a large enough number of species to conduct analysis). This allowed us to avoid the potential effect of small-ranged (or poorly sampled) species in determining the temporal dynamics of distribution and co-occurrence by, for example, increasing/decreasing species richness in particular time bins.
(d) Phylogenetic fields of species across time
Using the phylogenies and species' occurrences for each time bin, we estimated the phylogenetic structure of species co-occurrence within each focal range to determine the PFs of species. A species' PF represents the phylogenetic relatedness among species co-occurring within its geographical distribution, thus depicting the phylogenetic ‘neighbourhood’ represented by the community encountered by the focal species throughout its range . Therefore, PFs allow describing if species co-occur with closely or distantly related species or even with a random set of species throughout their ranges regarding the species pool present at each time bin. PFs were described with the PSV (phylogenetic species variability) index [16,44] that varies between 0 (reduced variability, clustering) and 1 (maximum variability, overdispersion), providing a single value of phylogenetic structure within species' ranges instead of an aggregate measure (e.g. mean) of phylogenetic diversity across sites.
For each time bin, individual PSV values as well as the average PSV of all species were computed. As before, we used these values to depict the temporal dynamics of PFs with box plots, considering all species in our database and also those species with more than seven occurrences. Non-parametric multiple comparisons were also conducted to statistically compare PF patterns among time bins. More specifically, to assess the temporal dynamics of PFs across the studied time interval, we focused exclusively on species present across all time bins (‘survivor’ species). For those species that were missing from a time bin but were otherwise present at the beginning and at the end of our studied time interval, we described their co-occurrences and associated phylogenetic structure by considering the intersection of the faunal lists of co-occurring species in the adjacent time bins where the focal species was present (i.e. time bins immediately before and after the missing time bin).
Observed PF patterns were compared with random expectations derived from 1000 runs of a null model. This null model simulated independent distributions among species in which a species' PF is generated by randomly sampling, without replacement, the observed number of co-occurring species from the time-bin phylogeny. The probability of sampling a co-occurring species was proportional to its range size (e.g. number of occupied cells, with large-ranged species being more likely to be sampled). This ‘range-based’ null model is appropriate for testing biogeographic patterns owing to the effect of species frequencies in co-occurrence metrics and phylogenetic structure [41,45] and the influence of large-ranged species in broad-scale biodiversity patterns [39,46].
(e) Macroevolutionary analysis
To associate PF patterns across time with their potential macroevolutionary drivers, we estimated speciation, extinction and diversification rates directly from the complete fossil occurrence dataset. We applied a recently developed Bayesian framework that estimates speciation and extinction rates and their temporal dynamics by using all fossil occurrences, including living species, and without relying only on first and last appearances [47,48]. This framework provides a robust probabilistic analysis that allows choosing the best fitting birth–death model under a Markov Chain Monte Carlo (MCMC) algorithm . Using the birth–death MCMC (BDMCMC) model allows modelling fossil occurrences as resulting from sampling (historical and geological conditions leading to the sampling, description and identification of fossils) and species diversification (speciation and extinction driving temporal changes in species richness). In addition, the BDMCMC model does not rely on arbitrarily predefined time bins and phylogenectic hypotheses, thus providing a continuous-time description of macroevolutionary rates explicitly considering extinct species (which is a desired property that contemporary phylogenetic methods based on molecular phylogenies generally miss ).
Using this framework, we estimated the number and magnitude of shifts in speciation, extinction and diversification rates of mammals across our studied time interval. Input data consist of fossil occurrences identified at the species level, with status (extinct or extant) as well as minimum and maximum ages  obtained from the fossil occurrence dataset. Our implementation of the BDMCMC model jointly estimated the number of rate shifts, their temporal placement, and the marginal speciation and extinction rates through time while considering the temporal heterogeneity of the preservation rate across species in our dataset. We also considered 10 randomly replicated datasets in order to accommodate uncertainty in the true ages of fossil occurrences . Finally, the estimated parameters were plotted in rate-through-time plots to visualize the temporal dynamics of the macroevolutionary process, combining the posterior estimates across all replicates. These analyses were conducted in the freely available computer program PyRate (; http://dsilvestro.github.io/PyRate/), using the default settings for the number of MCMC generations and sampling frequency. We excluded 40% of the samples as burn-in phase after checking for efficiency and mixing of the MCMC in Tracer .
(a) Species diversity and distribution across time
Species diversity and distribution varied noticeably across the six consecutive time bins analysed, resulting in considerable variation of species co-occurrence patterns. Considering our total 126 species with occurrence and phylogenetic data, species range sizes as well as total richness within time bins were higher in the more recent time bins when compared with the older (table 1 and figure 1a). However, instead of finding a concomitant pattern of higher species co-occurrence towards the Recent, species tended to co-occur with less species in the middle time bins (from 2 to 0.5 Ma) in comparison to the older (3–2.5 Ma) and youngest (0.5 Ma to the Recent) bins (table 1 and figure 1c).
At the level of species pairs, association between pairs was also lower in the middle time bins (figure 1e). Most of these patterns were even more apparent when considering only those species with more than seven occurrences in our database (figure 1b,d,f). Statistical comparisons among time bins showed significant differences for all these patterns (Kruskal–Wallis χ2 = 37.09, d.f. = 5, p < 0.001 for range size; Kruskal–Wallis χ2 = 42.39, d.f. = 5, p < 0.001 for within-range richness; and Kruskal–Wallis χ2 = 246.19, d.f. = 5, p < 0.001 for C-scores). Multiple comparisons among time bins showed that such differences were present especially between the youngest time bin and the middle ones but also between the older and middle time bins (figure 1; see the electronic supplementary material, tables S1–S3).
Temporal variation of species co-occurrences was also evident in the compositional structure of grid cells across time bins. The average proportional richness, which equals the reciprocal value of Whittaker's Beta diversity (1/Beta), showed a higher compositional turnover during the younger time bins (from 1.5 Ma to the Recent) with a lower number of species in an average grid cell and smaller average range size. Conversely, lower turnover and higher average richness and range size were observed towards the older time bins (from 1.5 to 3 Ma).
(b) Phylogenetic fields of species across time
Phylogenetic structure of species co-occurrence within individual species' ranges also varied across the studied time interval. PSV values revealed that average PFs of species within time bins varied slightly through time when considering all species (figure 1g), showing a clearer pattern when considering only those species with seven or more occurrences. More specifically, no significant differences among time bins was found in the former species set (figure 1g; Kruskal–Wallis χ2 = 6.7, d.f. = 5, p = 0.24), whereas significant differences were found for the latter species set (figure 1h; Kruskal–Wallis χ2 = 58.48, d.f. = 5, p < 0.001). Accordingly, lower PSV values were observed in the middle time bins when compared with the others (figure 1h). Thus, species seem to co-occur with more closely related species in those middle time bins when compared with co-occurrence with less closely related species within both the older and younger intervals. In general, across time bins, species with high within-range richness showed higher PSV values, whereas species with low within-range richness had lower PSV values. That is, species co-occurring with high numbers of species seemed to do so with more distantly related species, whereas species co-occurring with low numbers of species seemed to do so with more closely related species. However, observed patterns cannot be interpreted at face value and should be contrasted to expectations derived from the available species pool.
Comparisons between observed and null PF patterns showed that only a few species showed significantly clustered or overdispersed PFs (see the electronic supplementary material, table S1). Overall, across all time bins, most species with significant PFs showed an overdispersed co-occurrence pattern. Nevertheless, one generality derived from these comparisons was that most species with significant PFs appeared within the middle time bins, whereas most species showed random PFs at the oldest and youngest time bins (see the electronic supplementary material, table S4).
Focusing on the ‘survivor’ species, we found 20 species that were present across all time bins. These 20 species were analysed separately, providing a complete view of the temporal dynamics of species co-occurrences and phylogenetic structure. Similar to the overall patterns, most of these survivor species showed significant PF patterns within the middle time bins and random PFs at the older and younger time bins. Moreover, significant PFs were variable across species (table 2), showing highly idiosyncratic patterns without clear generalities regarding regions where species occurred or orders to which they belong (mostly associated with ecological specializations; see the electronic supplementary material, table S5). For instance, species from different regions (e.g. Europe and North America) as well as from different orders and having different ecologies (e.g. carnivores versus herbivores) showed both clustered and overdispersed PFs, specially within middle time bins.
(c) Macroevolutionary rates through time
The BDMCMC analysis supported a model with varying rates of speciation and extinction through the studied time interval. The model with the highest posterior probability included seven different rates of speciation and three different rates of extinction (see the electronic supplementary material, table S3). A model with constant rates of speciation and extinction was rejected, with posterior probability lower than 0.001. The rates-through-time (RTT) plots depicted the increase in speciation rate through time, with shifts at different points in time, specially those around 1.3, 1, 0.5 and 0.2 Ma with a higher than twofold increase from the older to the younger times and then decreasing towards the Recent (figure 2a). Conversely, the RTT plot depicted a highly constant extinction rate for most of the studied time interval, with shifts increasing this rate at around 1.2, 0.8 and 0.2 Ma when a 10-fold increase occurred towards the Recent (figure 2b). In balance, diversification rates also shifted, increasing around those same time points with a twofold increase around 0.5 Ma and then decreasing from around 0.3 Ma towards the Recent (figure 2c).
By integrating spatially explicit fossil data and phylogenetic information of species under a palaeomacroecological approach, we found that the combination of past climate changes and macroevolutionary processes determined the temporal dynamics of mammalian distributional and co-occurrence patterns. Across the time span considered here, we found higher significant variability of species co-occurrence and PF patterns within the middle time bins (that is some 2–1 Ma). Such patterns were roughly congruent with shifts in macroevolutionary processes and with known major climate changes around the same points in time. More interestingly, we also identified the idiosyncratic nature of species' responses to these drivers.
Palaeontological studies on the temporal dynamics of biodiversity have mainly focused on community-level patterns, describing changes in species richness and composition. Drivers of such patterns of community evolution—compositional changes due to originations, extinctions and migrations—are still hotly debated with some studies favouring climate changes (e.g. [51,52]) and others biotic interactions [53,54]. Alternatively, palaeontological studies focused on species-level patterns have concentrated on range shifts—changes in species' geographical distributions—associating them almost exclusively with climate changes [25,55,56]. In general, species are thought to respond in an individualistic manner to such changes and thus producing non-analogue communities . In agreement with this claim, our findings of species-level patterns of co-occurrence and phylogenetic structure support idiosyncratic responses of species to historical processes over time.
In terms of overall species co-occurrence, we found that, on average, total species richness within ranges was higher at the older and younger time bins, whereas lower richness within ranges was observed at the middle time bins. This pattern appeared in spite of an increase in species' range sizes and total species richness towards the Recent. In fact, such range size and richness increases towards the Recent may represent the outcome of uneven sampling , with younger time bins being better represented in the fossil record owing to their higher preservation potential . Accordingly, we could have expected a concomitant pattern of increasing species co-occurrence from the older to the younger time bins. However, this was not the case. We found similar co-occurrence patterns between older and younger time bins and lower co-occurrence patterns during the middle time bins. This lends conservative support to our findings over the potential effect of uneven sampling across time bins. Another potentially confounding factor could be the inclusion of large numbers of small-ranged (or poorly sampled) species in our analyses, since we considered all species with three or more occurrences in our database. However, analyses excluding those restricted species and focusing only on those with seven or more occurrences showed similar, even clearer, temporal dynamics in species distributional and co-occurrence patterns.
Species richness patterns are inherently linked to patterns in species range size, with both determining biodiversity patterns . Such relationships help explain our findings. For instance, lower species richness within ranges in the middle time bins (2–0.5 Ma) despite increasing overall richness in the system can be explained either by a reduction of species' range sizes and/or displacement to lower richness regions during those middle intervals. The first scenario, a reduction in species' range sizes, seems plausible given that species showed smaller range sizes in the older times when compared with the youngest time bin. However, lower species turnover along with a larger average proportional range size in the middle time bins runs counter to this scenario. Such lower species turnover within middle time bins is congruent with lower C-score values at these time bins, which implies higher aggregation between species pairs. Altogether, higher aggregation of species pairs and lower turnover can be further explained by the displacement of species, regardless of their range size, to low richness regions during those middle time bins.
Species shift their distributional ranges as a response to changes in climate . Such climate changes negatively affect the conditions where species occur, forcing them to disperse towards regions with more favourable conditions . Throughout our studied time interval, from the Late Pliocene to the Recent, there were several dramatic climate changes at the global scale that affected the biota, as have been shown for mammals . During the Plio-Pleistocene, strong mammalian faunal changes link to the beginning of major cooling and drying events happening around 2.5 Ma followed by a shift in the periodicity of climate cycles and the establishment of Pleistocene bipolar glaciations [28,59]. More specifically, important cold-shift oscillations particularly around 2 and 1 Ma have been associated with both community- and species-level changes such as great taxonomic and phylogenetic turnover as well as mixing of mammalian faunas and species shifting their ranges to track their original habitats (e.g. [15,26,55,56]). Indeed, it is around these time periods, which correspond to our middle time bins, where we found significant species-level changes in co-occurrence patterns such as the total within-range richness and species' aggregation patterns discussed above.
Changes in species' distributions inevitably result in changes at the level of local communities, especially in terms of species composition, thus producing non-analogue communities . Such compositional changes are further observed at the level of species by looking at the community formed by species co-occurring within individual ranges . Our PF approach allowed us to evaluate community evolution at the species-level by describing the phylogenetic structure of species co-occurrence within focal species' ranges . As previously stated, most studies of temporal biodiversity dynamics have focused at the community-level, mainly describing taxonomic turnover over time [28,60]. Only recently have these palaeontological studies begun to integrate phylogenies .
In a pioneering study, Raia  illustrated the applicability of the community phylogenetics approach  and its associated metrics to study the change in palaeocommunity phylogenetics. Raia showed, for Eurasian Ice-Age mammal communities, that temporal taxonomic replacement is not only due to species, but also clade turnover, with an overrepresentation of particular clades around 1 Ma resulting in phylogenetic clustering of those communities towards the Recent. Around that time, the full Ice Age was established and the so-called mammoth steppe had begun to dominate temperate environments, calling for habitat filtering forcing phylogenetically closed species to thrive in ‘cold’ habitats. In our case, when considering all of the species present at each time bin, we did not find differences in PF patterns among time bins. Alternatively, when excluding restricted species and focusing on those with seven or more occurrences, we found a tendency of species towards clustered PFs at the middle time bins, including the time around which the Ice Age began.
However, statistical comparisons showed that most species tended to present random PFs, whereas most species with significant PFs showed overdispersed PFs—co-occurring with distantly related species—around the middle time intervals (from 2 to 0.5 Ma) compared to the oldest (approx. 3–2.5 Ma) and youngest (0.5 Ma to the Recent) time bins. Therefore, contrary to Raia's findings , we found a tendency towards phylogenetic overdispersion of co-occurrence at the level of species ranges. Such discrepancy between Raia's  and our findings may represent their difference in scale, local and site-based versus regional and species-based, respectively. In turn, this would suggest either different processes or different responses to the same process at each scale, highlighting the need for site-based as well species-based analyses to understand the processes driving species co-occurrences [16,43,62].
Significant phylogenetic structure can result both from regional, historical events and local, ecological processes [61,63]. At broad spatial scales, species assemblages ultimately emerge from a combination of speciation, extinction and dispersal, namely historical processes [3,4]. However, ecological processes at local scales may still play a role in determining regional species assemblages [63,64]. For instance, Cardillo  showed that broad-scale assemblages of extant African carnivores are influenced by early rapid diversification followed by competitive sorting, with some assemblages being phylogenetically clustered and others overdispersed. In the same vein, Cooper et al.  found support for a tendency of extant mammalian assemblages to be phylogenetically overdispersed as a result of competitive exclusion among close relatives and/or ecological convergence among distant relatives. Indeed, determining the processes driving species assembly by using phylogenetic structure as a proxy is particularly challenging and may even be misleading . In our case, overdispersed PFs could be explained by competitive exclusion between closely related species and/or by evolutionary convergence of distantly related species, which would have in turn facilitated their co-occurrence [16,18]. A definite answer would require us to integrate data on traits mediating resource use . Nevertheless, our macroevolutionary analysis suggests that phylogenetic overdispersion of species co-occurrences may have resulted from increased speciation rates of different lineages at particular points in time, mainly around our middle time bins. In any case, our mixed findings of random, clustered and overdispersed PFs of species highlight their individualistic responses to climate changes and evolutionary processes assessed here, and perhaps also to ecological interactions, happening around those middle time bins.
Evaluating mammal species whose existence extends across the entire studied time interval further supports the notion for individualistic responses to the above-mentioned processes. Indeed, although these survivor species showed different PF patterns with overdispersed, clustered and even random fields, this fact highlights the idiosyncratic nature of species responses driving their co-occurrence patterns that, in turn, seems to be associated with relevant changes in climate and evolutionary processes taking place around the middle time bins. More interestingly, particular PF patterns (e.g. overdispersed fields) did not seem to be associated with specific characteristics of species. Carnivore species and herbivore species from different orders (e.g. Artodactyla, Perissodactyla and Proboscidea), for example, presented overdispersed, clustered or random PFs with some species, even changing from one pattern to another along consecutive time bins. Similarly, species restricted to particular regions showed all three PF patterns, except for North America, where species preferentially exhibited overdispersed patterns. This possibly depends on the nature of the landscape in North America, which is characterized by a large territory lacking internal natural barriers, which received intermittent influx of new species from the Old World through the Bering Strait during younger time bins. This, in turn, might have favoured strong competitive exclusion in North America as recently shown for carnivores .
Most changes in co-occurrence patterns of survivor species were observed during the middle time bins. All of these findings highlight the idiosyncratic nature of species responses and reveal a pattern of species co-occurrences similar to the evolution of local communities. That is, over time, species' PFs evolve as they encounter non-analogue assemblages throughout their ranges over time. Such non-analogue PFs (cf. ) may result from species' range shifts as well as from originations and extinctions continuously changing the species pool.
We showed here the application of a conceptual and methodological framework to investigate species-level patterns of co-occurrence and its phylogenetic structure using spatially and temporally explicit fossil data. We hope this framework helps to bridge the gap between neontology and palaeontology, offering a novel way to test old and new hypotheses on potential drivers of geographical biodiversity patterns through time. Although our findings should be interpreted at their considered phylogenetic, geographical and temporal scale [13,56,68], they highlight the importance of historical and evolutionary processes in shaping species co-occurrences. We provide here a first approximation to palaeomacroecological patterns at the species-level, hoping that future, more comprehensive analyses evaluating PFs through time will shed light on the ultimate determinants of the geographical patterns of biodiversity.
F.V. and J.A.F.D.-F. conceived the idea; F.V., F.C., P.R. and J.A.F.D.-F. designed the study; F.C. and P.R. compiled and curated the data; F.V. performed the analyses, coordinated the study, drafted the manuscript and led the writing. All authors contributed to the writing and gave final approval for publication.
We declare we have no competing interests.
F.V. was supported by a ‘Science without Borders’ BJT grant from CNPq. J.A.F.D.-F. is continuously supported by CNPq productivity grants. F.C. was supported by ‘STAR’ project grant.
We thank Tiago B. Quental, Thomas H. G. Ezard and Michael J. Benton for inviting us to contribute to this special issue. We are deeply thankful to Daniele Silvestro for thoughtful discussions and considerable advice on macroevolutionary analysis of fossil data. We also thank Matheus Lima-Ribeiro, W. Daniel Kissling, Tiago B. Quental and Sara Varela for discussions on the ideas behind this study. We thank Rebecca Crosthwait for kindly reviewing the English.
One contribution of 11 to a theme issue ‘The regulators of biodiversity in deep time’.
- Accepted December 9, 2015.
- © 2016 The Author(s)
Published by the Royal Society. All rights reserved.