Differential responses of Miocene rodent metacommunities to global climatic changes were mediated by environmental context

The study of how long-term changes affect metacommunities is a relevant topic, that involves the evaluation of connections among biological assemblages across different spatio-temporal scales, in order to fully understand links between global changes and macroevolutionary patterns. We applied multivariate statistical analyses and diversity tests using a large data matrix of rodent fossil sites in order to analyse long-term faunal changes. Late Miocene rodent faunas from southwestern Europe were classified into metacommunities, presumably sharing ecological affinities, which followed temporal and environmental non-random assembly and disassembly patterns. Metacommunity dynamics of these faunas were driven by environmental changes associated with temperature variability, but there was also some influence from the aridity shifts described for this region during the late Miocene. Additionally, while variations in the structure of rodent assemblages were directly influenced by global climatic changes in the southern province, the northern sites showed a pattern of climatic influence mediated by diversity-dependent processes.


Results
Metacommunities analyses. The cluster analysis identified four significant metacommunities (Fig. 1), which were identified with capital letters from A to D. Each metacommunity presented changes in the relevance of the different faunal components within its rodent associations (see Supplementary Table 2). Metacommunity A included fossil sites in which FC VI and FC V were dominant, sometimes accompanied by FC III. Sites included in metacommunity B were clearly dominated by taxa included in FC II. Metacommunity C was mainly dominated by FC III, with some relevance of FC II in the northern province. Finally, sites included in metacommunity D were restricted to the southern province in our database, and included higher percentages of species from the FC I and IV. The spatiotemporal dynamics of these rodent assemblages (Fig. 2) showed an initial dominance of metacommunity A in the Iberoccitanian Region with metacommunities B and C less represented. Progressively, during the Vallesian (around 11 to 9 Ma), the A metacommunity lost importance in favour of B and C in an incremental way. During the early Turolian (around 8.5 Ma), the metacommunity C was dominant in the Iberian Peninsula, and there was a northward shift of metacommunity B. Finally, the appearance of metacommunity D during the late Turolian (around 7.0 Ma) is associated to the north displacement of metacommunities B and C (Fig. 2). Geographical shifts or disappearance of these metacommunities in the northern province during the latest Turolian could not be assessed due to lack of record in this area between 6 and 5 Ma.
The nestedness analyses indicated a statistically significant nested pattern in the rodent assemblages from both bioprovinces, as well as in most of their different faunal components (Table 1). Diversity dynamics. The Shannon diversity index increased progressively in both provinces, although this increase was slightly more pronounced in the northern province (Fig. 3). Overall, this index exhibits fluctuations through time in both provinces (Fig. 3), but in an opposite pattern, with rough synchronicity between peaks in the northern province and valleys in the southern one.
According to the dynamic of rodent families, there were similar trends in both provinces, with early dominance of Cricetidae, Sciuridae and Gliridae (Fig. 4). Around 11-10 Ma the family Muridae is recorded in the Iberoccitanian region for the first time, and became the dominant family around 9 Ma in the southern province and 1 Ma later in the northern one (Fig. 4). The other families had much lower representation in the Iberoccitanian faunas, except for the families Castoridae, and Eomyidae in the northern province (Fig. 4). It is important to note the presence of the family Hystricidae only in the southern province, where it is recorded sporadically. Nestedness ranking of rodent assemblages was significantly correlated with taxonomical diversity only in the northern province (Table 2), which indicates an ordered pattern for variations of diversity in this province and suggests that the apparent loss/gain of genera in the northern province was diversity-dependent. The negative correlation between taxonomic diversity and nestedness ranking indicates that permanence of genera in face of environmental change within this province was higher if resources were shared among more equitably distributed genera (with similar proportions of species) rather than in communities with a higher incidence of closer relatives Environmental variables correlation. Our correlation analyses between climatic change (as measured by δ 18 O values) and the different variables studied here showed a high level of connection between environmental change and faunal dynamics. The Shannon diversity index was not correlated with δ 18 O in the southern province, although there was a significant correlation in the northern province, which shows that the general global cooling across the late Miocene appears to be related to increases in taxonomical diversity of its rodent faunas ( Table 2).
Proportion of species for most faunal components (FC) showed statistically significant correlations with isotopic values, which were slightly different in the southern and the northern province ( Table 3). Proportions of FC V and VI, both related to a set of middle Miocene genera ( Supplementary Fig. 2), significantly dropped in association with temperature increases in both provinces. FC III and FC IV, dominated by generalist species (Supplementary Fig. 2), correlated only marginally with isotopic variation in the southern province, while the former was significantly related to decreases in temperature in the northern province and the latter was recorded in only two northern fossil sites and statistical correlation was not possible. Finally, increases in the proportion of species of FC I, dominated by generalists and some species adapted to more open spaces ( Supplementary Fig. 2), were significantly related to increases in temperature in the southern province, but this was not the case in the northern province, probably because its absence in most of the northern rodent associations precluded a robust analysis. The nestedness ranking was significantly correlated with temperature in the southern province and non-significant in the northern province (Table 3). Regarding the independent FCs, in the southern province the FCs III, IV, V and VI were significantly correlated with the δ 18 O value, while in the northern province only the FC III was significantly correlated and correlations with FC V and FC VI were marginally significant ( Table 3). Looking into the sign of these correlations, the decrease of temperatures provoked gains of FC IV genera in the associations from the southern province, and gain of FC III genera in the northern province. In turn, the decrease of temperatures provoked a decrease in the fossil associations of genera included in the FCs III, V and VI in the southern province, and FC V and VI in the northern province (Table 3).

Discussion
Late Miocene rodent faunas from the Iberoccitanian region sharing similar ecological structures were classified into metacommunities, which followed non-random temporal and environmental assembly and disassembly patterns directly or indirectly related to environmental changes. Influence of climatic changes on the rodent communities proceeded through the integration of separate impacts on the different faunal components that integrated these assemblages. The proportions of species from each faunal component were significantly influenced by global temperature changes (Table 3). Additionally, our results strongly suggest that faunas from the northern and southern provinces showed differential responses to climatic changes; while temperature changes had a direct influence on the assembly-disassembly patterns of southern rodent faunas (Table 3), the influence of environmental changes appeared to be diversity-dependent in the northern province (Table 2); taxonomic diversity of rodent faunas in the northern province was influenced by temperature changes, while this was not the case for the southern province, where changes in humidity-aridity conditions might be more important 13,25 . This would be related to the effects of changes in temperature in the higher latitudes of the Iberoccitanian Region, where forested and more stable environments were dominant 7,36,37 . The temperate evergreen forests typical of the northern province in this period seem to be more affected by changes in temperature than in relative humidity 6,38 and, therefore, the rodent faunas adapted to these environments would be indirectly affected by landscape changes induced by temperature changes. Additionally, the increase in thermal seasonality associated with latitude increase made variations in temperature more relevant than relative humidity under these Miocene tropical conditions 6,39 , since it gave rise to the lack of specific food resources (e.g. fruits) during long periods of the year. The influence of diversity-dependent dynamics on the assembly-disassembly processes of rodent metacommunities in the northern province suggests that permanence of taxa in face of environmental change within this province was lower among close relatives than if they belonged to different genera. Although the functional role of diversity has been the subject of a long-standing debate in ecology, diversity tends to be correlated positively with ecosystem stability [40][41][42] , which is dependent on the differential response of species or functional groups to variable conditions, as well as the functional redundancy of species that have important stabilizing roles.
On the contrary, our results on the southern province suggest that, while the diversity of the local faunas was not correlated to changes in temperature (Table 2), the structure of the communities may have been influenced mostly by temperature changes ( Table 3). The significant correlations of the nestedness ranking with isotopic values showed that global temperature changes gave rise to an ordered substitution of the metacommunities and most of their components. However, some parts of that structure were possibly also determined by aridity changes. For example, there was no significant correlation between temperature changes and proportion of species in FC II (Table 3), which was dominated by forest adapted species (Supplementary Table 2). Due to the mostly dry climates of the southern province 19 , slight fluctuations of relative humidity had a crucial influence in forest environments and their associated fauna 9 .
In addition to the different ecological processes implied in the faunal dynamics of the provinces, faunal changes in relation to major environmental crises during the late Miocene appear to be delayed in the northern province when compared to the southern one (Figs 3 and 5). Eventually, changes in the ecological structure of rodent communities and turnover among metacommunities seem to be much more pronounced in the southern province than in the northern province (Fig. 5).
The simultaneous disappearance of metacommunities A and B and final dominance of metacommunity C in the southern province coincided with what has been called the Vallesian Crisis 43,44 , while there was a gradual substitution of metacommunities in the northern province during a much longer time interval during the late Miocene and started one million years earlier. Particularly interesting is the gradual change produced in the ecological structure of metacommunity B in the northern province across all the fossil sites studied (Fig. 5).
Casanovas-Vilar and Agusti 37 suggested that seasonality increase in the southern province was related to an exacerbated arid season while in the northern province there was simultaneous cooling and increase in aridity. These changes coincided with the gradual increase of thermal and hydric seasonality derived from the global cooling initiated in the middle Miocene 7,20,22,25,27,32 , which affected a set of species related to forested and humid  Additionally, the prevalence of metacommunity B in the southern province was probably related to humid pulses during the latest Aragonian (12-11 Ma) and early Vallesian (11-10 Ma) 9,18 , which allowed for the dispersal of forest elements from Central Europe through the northern province. The pulse of increment in taxonomical diversity around 10 Ma (Fig. 4) in the northern province is congruent with such a dispersal scenario and could be related to the first occurrence documented for murid taxa in this province 13 as well as the decrease in Gliridae (Fig. 4). While the initial concept of the Vallesian Crisis considered it as an abrupt event 44 and relevant changes have been observed in both macro-and micromammal faunas 7,20,32 , during the last years several authors have suggested that it could really be a succession of extinctions during a more prolonged time interval 7,24,27,32 . The coexistence of three different rodent metacommunities during the early Vallesian in the southern province as well as the prolonged persistence and gradual changes in the ecological structure of metacommunity B in the northern province appear to reinforce this interpretation of the Vallesian Crisis as a gradual event. Additionally, according to our results, the Vallesian Crisis rendered different outcomes in the southern and northern provinces of the Iberoccitanian Region, with the disappearance of metacommunities A and B in the former and the appearance of metacommunity C in the latter (Figs 2 and 5).
A second peak of diversity around 8.5 Ma in the northern province (Fig. 3) could be related to the displacement of metacommunity B towards higher latitudes within the northern province (Fig. 2), while metacommunity C became dominant in the Iberian Peninsula. The latter was mostly composed by murids, which were involved in a sudden diversification in Europe during the late Miocene 37,46 . This time was another moment of global cooling and aridification 13 , which favoured the survival and diversification of murids, probably due to their generalist condition. Nevertheless, regarding the southern province, taxonomic diversity in rodent faunas dropped around this date (Fig. 3), which was due to the relevance of rodent faunas highly dominated by numerous murid species (Fig. 4). The sudden reversion in the diversity trends in both the northern and southern provinces close to 8 Ma could be related to a recovery of previous climatic values (Fig. 4).
Noteworthy changes in rodent communities from the southern province around 7 Ma during the late Turolian were associated with the beginning of the Messinian, around one million years before the Messinian Salinity Crisis (Fig. 5). Metacommunity D appeared in the southern province (Figs 2 and 5) and completely displaced metacommunity C, which only remained in the northern province due to the arid conditions dominant in the   south and is involved in the last diversity decrease observed in the northern province (Fig. 3). Unfortunately, there is no record of rodent faunas in the northern province after this time period and the development of these metacommunities in that province is unknown at this moment. The sudden dominance of metacommunity D in the southern province could be related to the progression of decreasing global temperatures and aridity increase that, through glacio-eustatic sea-level changes, is at least partially responsible for the onset of the Messinian Salinity Crisis (MSC) around 6 Ma 22 . The MSC event represented a drastic increase in aridity, which affected the general pattern of faunal turnover in the Mediterranean region 7,22 and reduced taxonomical diversity of rodent faunas from the southern province (Fig. 3), largely dominated once more by murid species. Our work demonstrates that the influence of global and regional climatic changes on metacommunity dynamics is not simple and depends on the environmental characteristics of the affected areas as well as the characteristics of the species present in the community. In general, we found evidence of direct and indirect environmental controls on the variations of the ecological structure of rodent metacommunities, which were not only linked to temperature changes, but there was probably also an influence of the changes in aridity that occurred in this region during the late Miocene, particularly in the southern province. However, the diversity of rodent communities was more influenced by temperature than by humidity-aridity conditions in the northern province, through the effects of temperature and thermal seasonality changes on the dominant forest environment. Due to its buffer environmental conditions, the northern province of the Iberoccitanian region, or at least parts of it, played an important role as a humid fauna refuge where the metacommunities B and C survived successive aridity increases during the late Vallesian and the Turolian, before their alleged disappearance, which also evidences the stabilization effects of the higher diversity of the northern rodent faunas.
Although they were not drastic, we also found evidence of changes in the ecological structure of rodent faunas linked to the Vallesian and Messinian crises, particularly concerning the taxonomical diversity of the assemblages. We observed that diversity values were less variable in the southern province during the late Miocene, which is probably linked to the fact that aridity-prone faunas dominated rodent assemblages in this province. Notwithstanding, weak fluctuations in the Shannon diversity index through time in the southern sites were opposite to those in the northern province, which is consistent with the contrasting differences between these biogeographic provinces and their associated faunas.

Methods
Study area and fossil sites. The study area of this work spans the Iberoccitanian Region ( Supplementary   Fig. 1), comprising 117 fossil localities from the latest middle Miocene to the earliest Pliocene (12.6 to 4.9 Ma). This region is a key area area for the development of macroevolutionary studies due to the completeness, richness and abundance of fossil remains found in long-term and continuous stratigraphic sequences 14,47,48 . The study area exhibits strong environmental differences from the rest of Europe, which persisted in time due to its isolated position in the westernmost part of Europe 22,49 . Furthermore, there are two environmentally distinctive mammalian bioprovinces in this region 16,34 , recognizable since the Eocene 35,50 . The northern province includes fossil sites from the Rhône, Provence, Cucuron-Basse Durance and Languedoc-Rousillon basins from South-Eastern France, and the Vallés-Penedés basin from Catalonia (northeastern Spain). All the other fossil sites from the Iberian Peninsula are included in the southern province and are located at the Alfambra-Teruel, Alicante, Baixo Tejo, Castellón, Calatayud-Daroca, Duero, Fortuna, Granada, Guadix-Baza, Hijar, Murcia, Tajo, and Valencia basins (Supplementary Fig. 1).

Data base.
Our database is derived from Gómez Cano 6 and considers all rodent species recorded at the studied fossil sites, including 209 species. In order to reduce sampling biases, this database included well sampled fossil sites, which recorded a minimum of 100 m1 + m2 + M1 + M2, following van der Meulen and Daams 18 . Only a few fossil sites with a lower sample size were allowed as they were the only representatives of poorly sampled areas or because they were part of stratigraphically important sequences.
We employed this dataset of rodent species to compile a matrix with information on the number of species and percentage of species of each genus and family in each fossil site. As in other studies 16 , we only analysed taxa that were determined at the species level in each fossil site to avoid potential noise in the data due to unidentified taxa. We also assigned all the species in our matrix into the faunal components (FC) defined by Gómez Cano,et al. 7 , which grouped genera with similar ecological affinities ( Supplementary Fig. 2).
Finally, we also compiled independent bibliographic information for isotopic data contemporary to fossil sites included in our work 7 ; this work follows Gómez Cano, et al. 7 by fitting a smoothed curve to the isotopic information provided by Zachos, et al. 51 and interpolate δ 18 O values for the age of each fossil site. We used the global marine isotopic record because the continental isotopic record from the Iberoccitanian Region has low resolution in comparison to the rodent fossil record, although it is increasing at a steady pace [52][53][54][55] . Global isotopic record allow us to identify minimum ecological changes that affected rodent faunas.

Metacommunity identification.
In order to identify different metacommunities in the Iberoccitanian region during the studied time interval, we carried out a cluster analysis including the data from both bioprovinces.
Since we sought to group fossil sites according to the ecological affinities of their species instead of the taxonomical resemblance between them, we used faunal components as the studied variables 7 . These faunal components group together rodent genera with relatively similar ecological characteristics (as shown in previous works, see Supplementary Fig. 2 and references therein), which showed similar responses through time during the late Miocene 7 . Therefore, we calculated the percentage of species for every faunal component registered in each fossil site in relation to the total number of species in such site. The use of species percentages instead of the number of species allows the homogenization of all fossil sites for comparisons among them despite differences in species richness, which minimizes potential sampling biases (although they cannot be completely avoided).
The cluster analysis was calculated on Euclidean distances between group centroids and the clustering procedure was UPGMA (Unweighted Pair Group Method with Arithmetic Mean, also called Paired group) using the PAST 3.11 software 56 . Thereafter, we used the NbClust package of R software 57 to identify the significant number of groups within each cluster, which defined the number of metacommunities identified. Diversity dynamics. We calculated the Shannon index (H') to assess how diversity was related to changes in environment and ecological structure of rodent faunas within different metacommunities. This index takes into account the evenness of a dataset, so the cases with equitable numbers of entities between different types have higher diversity values. In ecology of modern ecosystems communities are customarily the studied cases, the types of interest are usually species and the entities of interest are commonly measured as number of individuals 58,59 . Nevertheless, since there are no available data on number of individuals for each species for all the studied fossil sites, and multiple taphonomic biases can affect to the representativity of proxies to individual abundance (such as number of identified specimens, inferred minimum number of individuals or number of molars), we used an approach that reflects taxonomical diversity rather than ecological diversity; we used the proportion of species (abundance of entities) for each genus (types) to calculate the Shannon index (H') for each one of the fossil sites (cases) in the two provinces. This measure gives an idea of evolutionary diversification within communities, establishing a gradient between communities in which resources are shared among a more or less equitable number of lineages. In order to study how this index changes through time we plotted the values for each one of the fossil sites against time and applied a local regression fitting (LOESS) over the data to visualize their trends throughout time. Finally, we also computed proportion of species in each faunal component as an additional biodiversity measure.
Since it is commonly proposed as a general model of community disassembly 60 , we also evaluated the possible presence of a nested structure in the assembly and disassembly patterns observed in the Iberoccitanian SCiEntiFiC RepoRts | (2018) 8:2502 | DOI:10.1038/s41598-018-20900-5 rodent faunas during the latest middle Miocene to Mio-Pliocene boundary. This model proposes that communities within disturbed systems exhibit nested structure such that the taxa included in poorer communities represent a confined subset of those in richer assemblages, rather than a random selection of those found in the entire species pool [61][62][63] . This pattern would imply that each taxon requires some minimal conditions to support population levels adequate to resist extinction, and that it can occur in all sites that attain these conditions. We calculated the nestedness of the rodent assemblages following the algorithm proposed by Rodríguez-Gironés and Santamaría 64 based on data of genus presence-absence matrices arranged by genus richness and number of occurrences. Following Furió, et al. 10 ,we used genera in order to avoid the possible noise derived from the species multiplicity within the fossil sites among sedimentary basins in both provinces 10 . Additionally, since species go extinct over a few m.y., in order to conduct the research over a long span of time, we have to use higher taxonomic units than species, and it is usually considered that congeneric species were ecologically similar. Finally, rodent taxonomy may change at the species level (e.g. difficulties in species differentiation due to anagenesis may become a problem in some lineages) but it is very consistent at the genus level.This algorithm calculates the nested subset temperature (a nestedness score) of each matrix in such a way that the lower the score, the more nested the structure of the community 65 . We calculated p values by means of a comparison to the distribution of scores generated by randomly shuffling the original matrices through 10000 Monte Carlo simulations (rows and sum totals were maintained constant). Nestedness analyses were run using the nestedness function as implemented in the R library Bipartite 57 and the null model 3 as suggested by Rodríguez-Gironés and Santamaría 64 , which is a constrained null model that accounts for the frequency of genera (column totals) and the genus richness of fossil sites (row totals) while sampling the null space uniformly, which minimises type I and II errors. Finally, when significant nestedness was identified, we compared the rank order in which assemblages were nested to their order based on diversity index and isotopic value using Spearman's rank correlation 66,67 . These analyses were performed for both bioprovinces as well as using independent matrices corresponding to the genera included in each faunal component 8 .

Influence of climatic change.
Finally, In order to evaluate the relevance of climate change in the establishment of the different metacommunities and their ecological and taxonomical characteristics (diversity and nestedness), we performed correlation analyses of the different measurements previously commented with the variations in the δ 18 O isotopic record as interpolated from Zachos, et al. 51 for the age of each fossil site included in this paper.