Body mass‐related changes in mammal community assembly patterns during the late Quaternary of North America

Pineda-Munoz, Silvia; Jukar, Advait M.; Toth, Anikó B.; Fraser, Danielle; Du, Andrew; Barr, W. Andrew; Amatangelo, Kathryn L.; Balk, Meghan A.; Behrensmeyer, Anna K.; Blois, Jessica; Davis, Matt; Eronen, Jussi T.; Gotelli, Nicholas J.; Looy, Cindy; Miller, Joshua H.; Shupinski, Alexandria B.; Soul, Laura C.; Villaseñor, Amelia; Wing, Scott; and Lyons, S. Kathleen, "Body mass-related changes in mammal community assembly patterns during the late Quaternary of North America" (2021). Faculty Publications in the Biological Sciences. 835. https://digitalcommons.unl.edu/bioscifacpub/835


Introduction
Classical niche theory suggests that trait-based approaches can be used to study the processes underlying patterns of association among species (Weiher andKeddy 1995, Kraft et al. 2008). Such processes include habitat filtering Jetz 2011), biotic interactions [e.g. competition (HilleRisLambers et al. 2012)] or dispersal limitation (Hubbell 2001). Under a habitat filtering model, species with similar habitat preferences or environmental tolerances are more likely to possess similar traits and co-occur because they have similar ecological requirements (Cornwell et al. 2006). In contrast, species with widely diverging habitat preferences or environmental tolerances co-occur less frequently and have different traits (Cornwell et al. 2006, Kraft andAckerly 2010). Biotic interactions, such as competition for the same resources, can cause species that have similar traits to exclude one another if competition is strong (Hardin 1960, Belmaker andJetz 2015). Additionally, coexistence of species may be constrained by dispersal limitation, such as geographic or climatic barriers (Condit et al. 2000, Hubbell 2001).
Previous studies have used null models (Harvey et al. 1983, Gotelli andGraves 1996) to study species co-occurrence patterns and, more recently, to classify the co-occurrence of each species pair as random, aggregated or segregated (Gotelli and Ulrich 2010, Blois et al. 2014, Lyons et al. 2016. Significantly aggregated pairs of species (aggregated pairs) will be found together across space more frequently than expected by chance, while significantly segregated pairs (segregated pairs) will co-occur less frequently than expected by chance. Pairs of species that are neither significantly aggregated nor segregated (random pairs) occur together with a frequency that is not different from that expected to occur by chance. Given the primacy of traits in mediating both environmental filtering and biotic interactions, functional traits can play an important role in community assembly and the type of association between a species pair (Ackerly and Cornwell 2007, Smith et al. 2016, Fitzgerald et al. 2017. For example, Smith et al. (2016) found that mammals within some trophic guilds, such as browsers, tend to form more aggregated pairs than mammals within other trophic guilds, such as carnivores.
Body mass is arguably one of the most ecologically informative functional traits of mammals because it is correlated with a number of traits such as diet (Andrews et al. 1979, Peters 1983, Pineda-Munoz et al. 2016. It is a relatively easy trait to measure in extant mammals and to estimate for extinct taxa. Moreover, there is a rich literature that uses body size to investigate mammal co-existence and community structure (Brown 1995, Bakker and Kelt 2000, Lyons and Smith 2013, Fraser and Lyons 2017. For example, Bowers and Brown (1982) showed that granivorous desert rodents of the same body size co-occur less frequently than expected by chance, an observation consistent with competitive exclusion. Similarly, Fraser and Lyons (2017) showed that mammal communities in northern temperate latitudes have higher dispersion in body mass compared to communities in the south, suggesting higher resource competition. Thus, body size as a functional trait offers an excellent means of studying the role of environmental filtering and biotic interactions, such as competition, in the assembly of mammal communities.

Climate change and biotic turnover in late Quaternary mammal fauna
The late Quaternary in North America is characterized by significant biotic turnover, as recorded by a well-documented mammal fossil record from 40 000 yr before present (ybp) to the present (Graham 1994, Graham et al. 1996. The Laurentide and Cordilleran ice sheets reached their maximum extent around 26 000 ybp, followed by a rapid deglaciation between 20 000 and 11 00 ybp (Clark et al. 2009, Denton et al. 2010. Species shifted their geographic distributions in response to these periods of climate change, sometimes resulting in non-analog communities (Graham et al. 1996, Lyons 2003, 2005. Humans also moved into North America during this time period, though the timing of their entry into North America is debated. Most research dates the dispersal of humans out of Beringia and into the rest of the Americas at ~13 000 ybp (Potter et al. 2018), where they altered ecosystems (Martin and Wright 1967, Alroy 2001, Faith and Surovell 2009). Both climate change and novel anthropogenic pressures on North American ecosystems have been implicated in the disproportionate extinction of mammal species weighing over 44 kg around 12 000 yr ago (Lyons et al. 2004, Koch and Barnosky 2006, Faith and Surovell 2009, Smith et al. 2018. Lyons et al. (2016) showed that both plant and animal co-occurrence patterns that were otherwise stable over the last 300 million years changed during the mid-Holocene, coinciding with increased human impacts. In particular, they observed a decline in the proportion of significantly aggregated species pairs that resulted in a lessening of the importance of biotic interactions in structuring co-occurrence of mammals on a continental scale (Tóth et al. 2019).

Co-occurrence patterns and body mass
Here, we use a trait-based approach to explore potential mechanistic reasons for the decline in aggregated species pairs from the Late Pleistocene to the present in North America (Lyons et al. 2016). In particular, we analyze the relationship between co-occurrence patterns and body mass difference (BMD) for pairs of species to understand the role of body mass in structuring species co-occurrence patterns. We use null models and niche theory to infer community assembly mechanisms and the reasons for the observed decline in aggregated species pairs. We investigate the following questions: 1) were certain pairwise body mass combinations more likely to change their association type over time, e.g. to change from aggregated to segregated? 2) Do different types of associations show different average BMD? 3) What is the relative importance of environmental filtering and biotic interactions (i.e. competitive exclusion) on mammal community assembly?
We address the third question by associating each combination of co-occurrence type (i.e. aggregations, segregations) and BMD (i.e. low, high) with a dominant ecological mechanism (Table 1). While not all species pairs are expected to be affected by the listed mechanism in Table 1, listed mechanisms are probably dominant if the group average BMD is lower/higher than random.
The re-organization of species distributions and community assembly patterns, the megafaunal extinction after the Late Pleistocene, and novel anthropogenic pressures make late Quaternary North American mammal communities an ideal study system for evaluating how contingencies can impact species co-occurrence patterns through time.

Dataset
We analyzed a dataset of pairwise associations produced in Lyons et al. (2016). In this work by Lyons and colleagues, species associations were evaluated using species-by-site presence-absence matrices from 369 localities across North America. The dataset covers three different time periods: Late Pleistocene (40 000 14 C-11 700 14 C ybp), Holocene (11 700 14 C-50 ybp) and Modern (50-0 ybp). Lyons et al. (2016) generated this co-occurrence data based on a collection of datasets known to include reliable data to minimize bias related to taxonomic resolution and taxa misidentification. The studied assemblages consisted of lists of species in a locality. Information about dataset collection and preparation for the analysis in Lyons et al. (2016) is fully detailed in Supplementary material Appendix 1.
For each species in Lyons et al. (2016), estimates of average body mass in grams were extracted from the MOM (mass of mammals) database (Smith et al. 2003) ver. 3. For extant species, estimates were averaged across sexes and species' geographic ranges. For extinct mammals, the MOM database compiled body mass estimates from the primary and secondary literature, and from regressions using tooth measurements. Body mass data were log10-transformed prior to analyses. BMDs were then calculated between each pair of North American mammal species. While some species' average body masses might have changed through time, our log-transformation makes it unlikely that results would be substantially affected; temporal changes in average body mass would need to change across an order of magnitude to impact our results. Lyons et al. (2016) performed the original evaluation of pairwise co-occurrence, whose results form the basis of our analysis. More details about their analysis can be found in Supplementary material Appendix 1, but the general logic of the method will be described here. Briefly, to determine whether each pair of species was aggregated, segregated or randomly associated for a given time period, a co-occurrence metric (i.e. C-score) was calculated for each pair of species in a species-by-site presence-absence matrix. A null distribution of C-scores was then generated by shuffling matrix elements, while preserving row and column totals (i.e. the 'fixed-fixed' algorithm). The observed C-score was then compared to the null distribution to determine whether a given species pair was significantly aggregated, segregated or random. These analyses were conducted using the software application PAIRS ver. 1.0 (Ulrich 2008). Overall, there were 25 459 species pairs across the three time intervals (Late Pleistocene, Holocene and Modern). Table 2 summarizes the results in Lyons et al. (2016). Previous research (Lyons et al. 2016, Tóth et al. 2019) has shown that these types of co-occurrence analyses are robust to differences in collection mode, temporal grain, spatial or temporal extent, taphonomic bias, taxonomic resolution and sampling biases.

Changes in proportion of association types through time
We counted the number of different types of pairs (aggregated, segregated and random) in each time interval (Table 2), and the proportion of pairs that changed from one association type to another (Fig. 1A). To evaluate the role of the extinct megafauna on the resulting patterns, we repeated the same calculations with species that were found in all three time intervals ('recurrent pairs'; n = 2370). These calculations Table 1. Inferred ecological mechanisms associated with each combination of pairwise co-occurrence type and body mass difference (BMD) level. Levels of BMD were determined by assessing whether the mean BMD for aggregated or segregated pairs was significantly lower or higher than that for random pairs. If aggregated or segregated species pairs have BMDs not significantly different from those of random pairs, this suggests neither environmental filtering nor limiting similarity were structuring species co-occurrence patterns regularly (at least with respect to body size). Such a finding may point to other mechanisms of importance (e.g. dispersal limitation, trophic interactions). Additionally, other indirect biotic interactions unrelated to BMDs may also play a role on community assembly.

Aggregations Segregations
Low BMD Environmental filtering: species occupying the same site and environment have similar body sizes due to similar ecological requirements and tolerances.
Limiting similarity: species with similar body sizes and resource requirements occupy different sites due to competitive exclusion. High BMD Limiting similarity: species occupying the same site cannot utilize the same resource and therefore have dissimilar body sizes.
Environmental filtering: species occupying different sites and environments have dissimilar body sizes due to dissimilar ecological requirements and tolerances.
* Table entries highlight ecological processes most relevant to this study, but other processes might be occurring.
excluded all mammal species (and associated species pairs found in Lyons et al. 2016) that went extinct from the Late Pleistocene to the Modern as well as rare modern species that are not found in the fossil localities included in the study.

Body mass differences (BMD)
To infer whether significant pairs were largely structured by environmental filtering or limiting similarity, we needed to determine whether the mean BMD of aggregated or segregated pairs was significantly larger or smaller than the mean BMD of the random pairs for each time period (Table 1). Because the number of species pairs within each association type differed by up to two orders of magnitude (Table 2; most species pairs were randomly associated), we created a null model to account for sample size differences among these categories. The null model was conducted for each pair of association type (i.e. aggregations, segregations, random pairs) within each time period. Focusing for now on a single time period and a single comparison between a pair of association types (e.g. aggregations versus random pairs in the Late Pleistocene), the null model was constructed as follows: 1) we calculated the absolute value of the difference in log 10 body mass (BMD) for each species pair in each association type. 2) We calculated each association type's mean BMD.
The difference between those means is the test statistic in our Step 3 was repeated 10 000 times to get a null distribution of differences in mean BMD between association types. 5) To compute a two-tailed p-value, we calculated the proportion of all the values equal to or more extreme than the positive and negative value of the observed difference in mean BMD (from step 2). All analyses were performed using R ver. 3.6.1 (R Core Team).
We did not specifically test for phylogenetic autocorrelation in the body mass differences analyses for two reasons. Firstly, the higher-level taxonomic composition of the dataset changes very little across the studied intervals (i.e. where we lose higher level taxa such Proboscideans, they do not make up the majority of lost diversity because they are not speciesrich in Late Pleistocene North America). Thus, we expect that the degree of phylogenetic signal for body mass is similar across all the time bins. Additionally, we get similar results when we repeat the analysis including only recurrent species pairs, which completely removes variability in phylogenetic composition across time intervals. Secondly, it is unlikely that our statistical comparison is phylogenetically autocorrelated. We compared the difference in mean BMD between association types to a null distribution. The data points are not individual species and we compare a single value (mean BMD) to a distribution of randomized values. Additionally, by taking the difference, we remove some phylogenetic nonindependence. Though the formation of the underlying species pairs may reflect phylogenetic signal of body mass and habitat preference (and, probably, a correlation between the two), the means we compare are unlikely to be phylogenetically autocorrelated (Felsenstein 1985, 1988, Garland Jr et al. 1992, Grafen 1992.

Changes in proportion of association types through time
Based on the PAIRS analysis of Lyons et al. (2016), the absolute number of pairs for each association type differed by up to two orders of magnitude (Table 2), as most species were randomly associated in all time intervals (Table 2). Across all time intervals, North American mammal communities experienced a decline in the proportion of significantly aggregated associations and an increase in segregated associations relative to the total number of significant pairs from the Late Pleistocene to the Modern. This is primarily due to random pairs becoming segregated (Fig. 1). In the Late Pleistocene, 59% of significant pairs were aggregated. Aggregations fell to 46% of significant pairs in the Holocene and then to 41% in the Modern. Many species pairs changed, becoming, for example, segregated when they had previously shown no significant association or disassociation (Fig. 1). However, it was rare for aggregations to become segregations and vice versa; most pairs tended to transition from aggregated to random or from random to segregated from the Late Pleistocene to the Holocene (Fig. 1). We observe the same pattern when we look only at recurrent species pairs (Table 2, Fig. 1B).

Body mass difference across time intervals
The results of the null model show that, in the Late Pleistocene, the average BMD for aggregated species pairs was significantly smaller than for segregated pairs (p < 0.05) and random pairs (p < 0.05). In the Holocene, the BMD of segregated pairs declined, becoming significantly lower than for random pairs (p < 0.05). Differences in mean BMD between random and aggregated remain significantly different. In the Modern, neither aggregated pairs nor segregated pairs have significantly different BMD from random pairs.
During the Late Pleistocene and Holocene, the average BMDs of aggregated and segregated species pairs that are recurrent in all time intervals were significantly smaller than those of random pairs (p < 0.05). In the Modern, neither aggregated pairs nor segregated pairs have significantly different BMD from random pairs (Fig. 2, Supplementary material Appendix 1 Fig. A2).

Body mass and association type
Certain body mass combinations are more likely to change association type than others (Fig. 3). Randomly associated pairs of species within the small body mass categories (10 2 and 10 3 g) increased their proportion in the Holocene compared to the Late Pleistocene (dark-grey cells in Fig. 3A). Meanwhile, the number of random species pairs that included larger-bodied mammals (in the 10 5 and 10 6 g categories) decreased during the same time interval (lighter gray cells in Fig. 3A, Late Pleistocene versus Holocene). We see no change in the relative proportion of random species pairs per body mass combination in the analysis including only recurrent species pairs.
From the Late Pleistocene to the Holocene aggregated species pairs with body masses of 10 2 -10 3 g increase and species in pairs with larger body masses (10 4 -10 6 g; dark-blue Fig. 3A) decrease. Segregated pairs show an increased proportion of species in the medium size range (10 3 and 10 4 g, dark-red cells Fig. 3A). Bigger mammals (10 5 g) experience a decrease in the segregated pairs. Across the Holocene-Modern transition, random pairs show little change in proportional abundance of body sizes (Fig. 3A). This is true using all species pairs or when restricted to recurrent pairs. Aggregated species pairs with body sizes of 10 4 -10 6 g slightly increase in relative proportion in the analysis with all species pairs from the Holocene to the Modern (blue cells Fig. 3A). Relative abundance of smaller species in pairs (10 2 and 10 3 g) decrease. Pairs with species of a body mass of 10 2 -10 4 g formed an increasing number of segregations from the Holocene to the Modern in both the analysis with all species pairs as well as with only recurrent species pairs (red cells Fig. 3A and B, Holocene versus Modern).

Discussion
In the Late Pleistocene, significantly aggregated pairs of species had lower average body mass differences (BMD) than random pairs (Fig. 2). Similarly, the average BMD of segregated pairs was higher, though not significantly, than the average BMD of random pairs (Supplementary material Appendix 1 Fig. A1). Environmental filtering therefore appears to have exerted an important control on species co-occurrence during the Late Pleistocene (Table 1): species with similar body mass values likely aggregated as a result of similarities in environmental tolerance and habitat preference (Danell et al. 2006, Bakker et al. 2016. For example, the bighorn sheep Ovis canadensis and the mountain goat Oreamnos americanus were significantly aggregated in the Late Pleistocene. These two species have similar body sizes [74 644 g and 72 500 g respectively; Smith et al. (2003)] and occur at cold, highelevation habitats (Nowak 1999). Previous research has suggested that resource partitioning might allow bighorn sheep and mountain goats to occur in the same habitats and avoid competition, despite ecological similarity (Reed 2001). Thus, high resource availability may have allowed many species that share climate tolerances to coexist during the Late Pleistocene.
In the Holocene, we observe a decline in percentage of aggregations and increase in percentage of segregations, the latter of which was concentrated among the medium-sized mammals (10 2 -10 4 g Fig. 3). The Pleistocene-Holocene transition was characterized by range shifts and expansions amongst a majority of mammals (Lyons 2003, Tóth et al. 2019, resulting from weakened climate gradients Figure 2. Distribution of species' pairwise body mass difference (BMD) for each association type for the Late Pleistocene, Holocene and Modern intervals for (A) all species pairs and for (B) species pairs that can be found in all time intervals (recurrent species). p-values were calculated using a null model analyzing body mass differences (BMD) (see Body mass difference model in Methods and Supplementary material Appendix 1 Fig. A1, A2). Seg: segregated; Rand: random; Agg: aggregated. BMD (body mass difference) = |log 10 (body mass (g) sp.1) -log 10 (body mass (g) sp.2)|. Figure 3. Number of species pairs relative to the total number of pairs per association type and time interval in each pairwise combination of body mass log 10 order of magnitude (g) for (A) all species pairs and for (B) species pairs that can be found in all time intervals (recurrent species) as extracted from Lyons et al. (2016). Darker colors indicate higher proportions as shown in the legend. N is the number of species pairs in each time interval and association type. BM: body mass. (Dean et al. 1984, Shuman andMarsicek 2016) and landscape-scale changes in the wake of the megafauna extinction (Malhi et al. 2016). Holocene range shifts likely had complex effects on the numbers of segregated species pairs. These segregations could result either from range contraction or decreased occupancy. However, this is not borne out in previous analyses, which show increased occupancy and range sizes in the Holocene (Tóth et al. 2019). On the other hand, the increase in segregations during the Holocene may result from the loss of biotic associations and enhanced importance of abiotic segregations such as habitat preferences during the Pleistocene-Holocene transition (Tóth et al. 2019).
Changes in small mammal communities may also explain the increasing number of segregated pairs. Although Late Quaternary extinction events mostly affected larger mammals (species in the 10 5 -10 7 g body mass categories in Fig. 3A), communities of smaller mammals also experienced a marked decline in richness and evenness in the Holocene (Grayson 2000, Blois et al. 2010, Lyman 2014. While human impacts such as landscape modification and overhunting were likely the major driver for the Late Pleistocene-Holocene megafaunal extinction (Martin and Wright 1967, Alroy 2001, Lyons et al. 2004, Faith and Surovell 2009, Smith et al. 2018, Stephens et al. 2019, the changes in North American small mammal communities began prior to first human settlements, during the Late Pleistocene (Grayson 2000, Blois et al. 2010, Lyman 2014, implicating climate change (Dean et al. 1984, Shuman and Marsicek 2016) in their observed decline in richness and evenness. Generalist species that thrive in disturbed environments such as Peromyscus spp. doubled their abundances, and perhaps outcompeted more specialized small mammals during the Holocene (Blois et al. 2010).
Geographic dispersal limitations can also play a key role in the observed species' co-occurrence patterns (Hubbell 2001, Allen et al. 2006. Barriers to dispersal such as habitat fragmentation and land-use change, which progressively increased during the Holocene (Stephens et al. 2019), may prevent some species from occurring in sections of their former ranges and cause the observed increase in segregations and decrease in aggregations in small mammals (Ordonez et al. 2014, Haddad et al. 2015. Other studies have shown that habitat fragmentation can prevent the co-occurrence of species that normally coexist, and that this especially affects small mammals (Cáceres et al. 2010).
Despite changes in the number of significant pairs during the Holocene, the average BMD of aggregated pairs was significantly smaller than in random pairs, as in the Late Pleistocene. This indicates that similar ecological processes, including environmental filtering continued to structure species aggregations in Holocene ecosystems (Table 1). Segregated pairs, however, had a significantly lower average BMD relative to random pairs. Specifically, those species pairs that became segregated in the Holocene but were randomly associated in the Late Pleistocene mostly belong to species pairs in the 10 2 -10 4 g body mass categories (dark red cells Fig. 3A-B). This limited body mass similarity among co-occurring species suggests that competition played a larger role in the Holocene than in the late Pleistocene. Increased segregations among species in these body size categories could be explained by competition that limited body mass similarity (Table 1), together with limiting similarity in other traits that correlate with body size such as diet preference (Pineda-Munoz et al. 2016). Previous studies have shown that species of similar body mass rarely co-exist when habitat disturbance is high (Bowers and Brown 1982, Brown and Bowers 1985, Brose 2010, Belmaker and Jetz 2011, 2015 and the Late Pleistocene-Holocene transition was a time of rapidly oscillating climate and habitat change (Dean et al. 1984, Shuman andMarsicek 2016). Thus, climate fluctuations could have increased the role of ecological processes such as competitive exclusion, which is represented by the formation of segregated species pairs.
When we remove the large-bodied extinct species from the analysis, however, both Late Pleistocene and Holocene segregated pairs have significantly lower BMD than random pairs (Fig. 2B), indicating that the same body mass-related ecological processes (e.g. competition-driven segregation of ecologically similar species) drove species segregations among surviving species in the Late Pleistocene and Holocene (Table  1). Formation of aggregated and segregated pairs among surviving mammals during the Late Pleistocene and Holocene appear to have been driven by the same ecological processes of environmental filtering and competition, respectively. Though range expansions and shifts were widespread among Holocene mammals (Lyons 2003, Tóth et al. 2019), they did not lead to changes in the ecological rules governing species pair formation.
During the modern interval, however, the mediating effect of body mass as a driver of community assembly processes such as environmental filtering and competition appears to have been diminished relative to the Late Pleistocene and Holocene. Although a greater proportion of pairs overall are significant after the Holocene-Modern transition (Table 2, Fig. 3A and B), we observe no significant differences in average BMD across association types (Fig. 2). We suggest that the increased number of both aggregated and segregated pairs as well as decreased importance of body mass in community assembly relate to the acceleration of anthropogenic impacts during the modern interval.
The 20th century saw the fastest and most spatially extensive human landscape modification, including the transformation of large portions of the landscape into farmland, expansion of human transportation networks, and rapid growth of urban environments (Ellis et al. 2013, Boivin et al. 2016. The result has been significant fragmentation of natural ecosystems (Cáceres et al. 2010, Haddad et al. 2015. The proportion of total pairs that were segregated increased from 4.35% to 7.34% from the Holocene to the Modern, likely reflecting enhanced geographic range fragmentation and reduced range overlap among small mammals in particular (Table 2) (Terry andRowe 2015, Terry 2018). In contrast, we observe an increase in aggregated pairs among larger species (10 3 -10 5 g) (Fig. 3). Large-bodied species are less dispersal limited than smaller-bodied species (Brown 1995, Lyons 2005 and less affected by habitat fragmentation resulting from the expansion of human habitation and transportation networks, as well as the patchy distribution of protected lands (Hubbell 2001, Allen et al. 2006. Future directions should include exploring the role of habitat variation across space on driving disjunctive ranges and structuring species co-occurrence patterns across the Holocene. Humans also transplant species on which they rely for food and sport, effectively expanding their geographic ranges, including the coyote, whose range increased many times over due to human sport hunting preferences (Hill et al. 1987). The cumulative effects of agriculture and species transplantation were, therefore, increases in range size (Lyons 2003, Lyons et al. 2010, Tóth et al. 2019) that may have resulted in enhanced numbers of species aggregations among large mammals (Fig. 3). The observed increase in the relative proportion of aggregated pairs could relate to increased range overlap due to geographic range expansion (Lyons 2003, Cáceres et al. 2010, Lyons et al. 2010, which may cause species to cooccur more often than they did in the past. Beyond being a pattern we observe, range expansion may also cause biotic homogenization, an ecological process by which ecosystems lose their taxonomic, genetic or functional uniqueness (Olden andRooney 2006, Davey et al. 2012). The ecological consequences of human landscape modification are therefore likely to have been far reaching.

Conclusions
We find that the role of body size in mediating species associations changed over the last 40 000 yr. During the Late Pleistocene and Holocene, aggregated and segregated species pairs were more similar in body mass than random pairs (Fig. 2). In the Modern, however, differences in body size were indistinguishable from random pairs for aggregated and segregated pairs, suggesting a reduced importance of the ecological processes that drove pair formation in previous time periods (Table 1), as can be inferred from body size.
Our results show that average BMD between species pairs can be used to understand community assembly patterns and mechanisms in modern and fossil ecosystems. During the Late Pleistocene of North America, aggregated species pairs had smaller average BMD compared with random pairs, which suggests environmental filtering was an important process for structuring communities. We observe the same pattern during the Holocene in significantly aggregated pairs of species. Segregated pairs, however, have significantly lower BMD in the Holocene compared with segregated pairs in the Late Pleistocene. While environmental filtering was still important, competition may have played an increased role in structuring Holocene communities. In modern mammal communities, however, we observe no significant effects of body size on species co-occurrence from comparing the BMDs of significant pairs to the BMDs of random pairs. Our analyses suggest that human impacts may have fundamentally altered the way mammals coexist on the North American landscape.

Data availability statement
Data available from the Dryad Digital Repository: <https:// doi.org/10.5061/dryad.pg4f4qrmw> (Pineda-Munoz et al. 2020). Data can be accessed from the Supplementary material. Data includes all analyzed pairs of species; whether they occur in the Late Pleistocene, Holocene or Modern; type of association; and body mass of each species in a pair. We also include the R code for the null model analysis.