Neogene precipitation, vegetation, and elevation history of the Central Andean Plateau – Science Advances

Abstract

Andean uplift played a fundamental role in shaping South American climate and species distribution, but the relationship between the rise of the Andes, plant composition, and local climatic evolution is poorly known. We investigated the fossil record (pollen, leaves, and wood) from the Neogene of the Central Andean Plateau and documented the earliest evidence of a puna-like ecosystem in the Pliocene and a montane ecosystem without modern analogs in the Miocene. In contrast to regional climate model simulations, our climate inferences based on fossil data suggest wetter than modern precipitation conditions during the Pliocene, when the area was near modern elevations, and even wetter conditions during the Miocene, when the cordillera was around ~1700 meters above sea level. Our empirical data highlight the importance of the plant fossil record in studying past, present, and future climates and underscore the dynamic nature of high elevation ecosystems.

A major uncertainty in climate modeling is how future precipitation will change over the Central Andean Plateau (CAP) and Amazon Basin (1). While some model outputs suggest feedbacks that will lead to massive drying (2, 3), others suggest wetter-than-modern conditions (46). Ongoing changes in temperature and atmospheric CO2 concentration (pCO2) will, by the end of this century, provide similar conditions to those found in the late Miocene and early Pliocene (7). Between ca. 10 million years ago (Ma) (late Miocene) and 5 Ma (early Pliocene), the northernmost part of the CAP (NCAP) doubled in elevation (8, 9). Mechanistic models (e.g., RegCM3 and ECHAM) project that this uplift increased precipitation over the CAP, creating essentially modern climate patterns (Fig. 1) (5, 911). The CAP paleobotanical record provides physical evidence that can be used to reconstruct paleoprecipitation, providing elements for testing the validity of climate models of the geologic past and for refining predictions of future climates. Here, we provide both Miocene and Pliocene paleobotanical evidence to document the effects that uplift of the NCAP had on regional ecosystems and climates.

Geographic location of the CAP, summary of paleoprecipitation estimates from previous studies. (A) Topographic map showing the location of the CAP and the location of the Descanso-Yauri Basin (black triangle) in the Peruvian Andes. Gray triangles indicate the geographic location of sites from the modern palynological dataset and their environmental distribution (elevation and mean annual precipitation). m.a.s.l., meters above sea level. (B) Simulated precipitation based on the ECHAM model when the CAP had 100% (600 to 700 mm/year) and 40% (100 to 200 mm/year) of its current elevation [figure taken from (9)].

The modern NCAP sits at ~4000 m of elevation and has a mean annual temperature (MAT) ~8 to 9C and a mean annual precipitation (MAP) ~500 to 800 mm. The CAP is characterized by cold and strong winds throughout the year, and it experiences extreme diurnal and nocturnal temperature changes reaching up to 30C (12). Most precipitation in the CAP falls during the austral summer as a result of easterly winds bringing moisture from the Atlantic Ocean via Amazonia (13). The modern high-elevation vegetation types of the NCAP are grasslands and shrublands, known locally as puna.

Here, we studied the Neogene paleobotanical record from the NCAP (Fig. 1) aiming to (i) estimate regional paleoelevation and paleoprecipitation and (ii) reconstruct its paleovegetation. The plant fossil material was collected in the Descanso Formation, Descanso-Yauri Basin, southern Peru. This basin has an approximate area of 2000 km2. We studied two members (or lithostratigraphic units) from the Descanso Formation that were previously dated using radiometric methods: Member B, dated to ca. 18.7 to 9.1 Ma (middle to late Miocene) (8, 14), and Member C, dated to ca. 4.8 to 3.9 Ma [Fig. 2; early Pliocene; (8)]. Exposed along most of the basin, Member B is the thickest unit of the Descanso Formation and was deposited in a braided-river environment (8). There is an angular unconformity (time gap) between members B and C of approximately 4 Ma during which an uplift of ca. 2500 m has been inferred on the basis of isotopic data (8). Member C is a thin unit deposited in a fluvial-lacustrine system at a similar elevation as it is found today (8, 15).

(A) Geologic map for the Descanso-Yauri Basin, showing where the samples were collected [modified from (8)]. (B) Stratigraphic column for El Descanso Formation, including members A, B, and C [modified from (8)]. Bars next to the samples correspond to the uncertainty in the stratigraphic depth for the sample or group of samples. Undulating line represents the unconformity between members B and C. (C) Photographs of some of the macro- and microfossils collected for this study. From top to bottom: Fossil log with 7 m of length, legume fossil wood anatomy (scale bar, 1 mm), Ribes sp. leaf (scale bar, 10 mm), and Podocarpus sp. pollen (scale bar, 10 m). Photo credits: First and second rows: C. Martnez, Cornell University and Smithsonian Tropical Research Institute (STRI); third row: A. Aliaga, Museo de Historia Natural de Lima; fourth row: J. E. Moreno, STRI.

The plant fossil record from the Descanso-Yauri Basin studied here included palynological (pollen and spores) and macrofossil samples (permineralized wood and compressions and impressions of leaves and fruits) of Neogene age (Fig. 2). The age of the fossiliferous localities was determined using previously published stratigraphic information (8), together with new zircon U-Pb geochronologic analyses of two tuff samples from the Member B that yielded ages of 10.03 0.16/0.18 Ma and 12.07 0.73/0.74 Ma (fig. S1 and the Supplementary Materials). A total of 187 paleobotanical samples were collected from 88 localities from members B and C. Seventy-five palynological samples were analyzed, containing 168 palynomorph taxa in 5389 individual occurrences (table S1 and the Supplementary Materials for raw counts). From those, 154 palynomorphs were described (see the Supplementary Materials for descriptions and plates for each palynomorph). Thirteen samples of permineralized wood were assigned to four morphotypes. Ninety-nine samples consisted of macrofossil compressions and impressions of leaves and fruits and were assigned to eight morphotypes (see the Supplementary Materials for descriptions and plates for each macrofossil morphotype).

Wood samples were preserved as silica permineralizations and were only found at the top of Member B (Fig. 2). Eight specimens had taxonomic affinities with Fabaceae and one with the fossil genus Anacardioxylon (Anacardiaceae). Within Fabaceae, two samples were identified within the Tribe Ingeae, five samples within the mimosoid clade, and one as the fossil genus Andiroxylon sp. (see the Supplementary Materials). Tree height was estimated for the Andiroxylon sp. sample based on its diameter, and estimates ranged from 29.8 to 34.6 m. In tropical high-elevation ecosystems, plant growth is dependent on the number of hours of freezing per day (16). A correlation between tree height and elevation across six modern Andean elevation transects shows that trees taller than 30 m are not found above ~1100 m of elevation in the Andean region (17). Nevertheless, assuming that plants shift their distribution upward when the temperature increases (18), a temperature correction for the late Miocene would produce higher estimates. On the basis of these observations, we infer that trees of the size found in Member B would not be found above ~1750 m during the late Miocene. This provides a maximum likely paleoelevation for Member B at the time of deposition.

The traits of woody plants found in Member B are consistent with a warm and wet setting. All wood samples had simple perforation plates. Twelve of the 13 samples collected had diffuse porous wood and nondistinct growth rings. Six of the seven samples identified had large vessel diameters (>150 and up to 270 m) and high proportions of axial parenchyma cells. Diffuse porosity, nondistinct growth rings, vessels of large diameter, simple perforation plates, and a large proportion of axial parenchyma cells were dominant wood anatomical characters in our fossil assemblage. That combination of characters is commonly found in trees of large size, having lowland affinities, and growing in sites with high MAP regimes and low seasonality (19, 20).

Leaves were found in both members; Member B only had palm leaf fragments (Arecaceae), whereas Member C had leaves with taxonomic affinities to five genera and one family that are today present in the Puna: Ribes (Grossulariaceae), Berberis (Berberidaceae), Polylepis (Rosaceae), Polystichum (Dryopteridaceae), Equisetum (Equisetaceae), and Juncaceae (see the Supplementary Materials for descriptions and photographs of each morphotype). The palynological record shows that the floristic composition of Member B is dominated (>50 counts) by the families: Podocarpaceae, Cyatheaceae, Polygonaceae, Poaceae, Polypodiaceae, Chenopodiaceae, Chloranthaceae, Malvaceae, Asteraceae, Lycopodiceae, Araceae, Solanaceae and Caryophyllaceae, while Member C is dominated by Poaceae and Cyatheaceae (table S3). Studies across modern Andean transects have shown that altitudinal vegetation gradients are accurately reflected by modern pollen rain data and thus that fossil palynological assemblages can be used to reconstruct aspects of past Andean ecosystems (21).

Taking the sum of the paleobotanical data, we developed a new quantitative coexistence analysis to reconstruct paleoprecipitation and paleoelevation. This coexistence analysis uses climate and elevation information associated with the nearest living relatives of the identified palynomorphs and macrofossil taxa to estimate a mutual climate and altitudinal range (table S3). Although climate tolerances of species could have shifted since the Miocene, no substantial changes seem likely to have occurred given physiological uniformitarianism (22). Two distribution datasets (palynological and macrofossil) were described through bivariate probability density distributions of modern taxon occurrence along elevation and precipitation gradients. Global temperatures during both, middle to late Miocene and early Pliocene times, were warmer than preindustrial temperatures. Therefore, we performed a correction of temperature estimates that resulted in average displacements of the elevation by ~650 and 440 m for members B and C, respectively, using modern lapse rate for the Central Andes (fig. S3) (23). Each distribution dataset (palynological and macrofossil) was first analyzed independently, and later after a mixture sensitivity analysis was applied (fig. S4), final estimations were weighted equally for precipitation and for the elevation of Member C (1 = 2 = 0.5). However, an unequal weight was assigned for elevation estimates for Member B due to the high sensitivity found in our analysis depending on the dataset used (1 = 0.615). For Member B, the median elevation based on the mixture model (palynological and macrofossil data) was 1636 m (interquartile range, 971 to 2850 m), whereas for Member C, the median estimated elevation was 3780 m (interquartile range, 3340 to 4140 m). These corrected elevation estimates are congruent with previous paleoelevation estimates based on isotopic data for the NCAP (fig. S5) (8). Median annual precipitation for Member B based on the mixture model was 1671 mm/year. (interquartile range, 1265 to 2077 mm/year), whereas for Member C, it was 862 mm/year (interquartile range, 620 to 1121 mm/year; Fig. 3). These estimates of environmental space show a trend that progressed in time toward ascending elevation and decreasing precipitation (Fig. 3). An additional comparison with environmental space from Holocene records from the CAP shows that changes in the environmental space from Miocene to Pliocene were much greater than those transitioning from Pliocene to Holocene (Fig. 3). It is worth noting that nonanalog assemblages can still be used to provide quantitative paleoclimate reconstructions. Analogs rely on the co-occurrence of species in their realized niche space, but the fundamental niche (sensu Hutchinson, 1957) (24) of species is broader than their observed realized niche. By using the range tolerance of the genera, probabilities can be generated for overlapping environmental characteristics, even if the genera do not co-occur today (25).

Estimated environmental space for the middle to late Miocene (Member B; gray), the early Pliocene (Member C; red), and the Holocene (blue) at Descanso-Yauri Basin, showing probability density functions for precipitation (right) and elevation (bottom). These estimations are based on a mixture of micro- and macrofossil sets.

Our reconstructions for the middle to late Miocene indicate that when the southern Peruvian Andes were at half of their modern elevation (~1600 m), MAP was about three times higher (~1700 mm/year) than today. In contrast, regional climate models (ECHAM, GENESIS, PLASIM, and RegCM3) predict that half-modern CAP elevations would lead to much-drier-than-modern conditions (3, 5, 10, 11, 26, 27). For example, the global ECHAM model indicates that when the Andes were at 40% of their modern elevation, MAP was 100 to 200 mm, or about ~80% drier than today (Fig. 1B) (9). Previous paleobotanical studies also suggest that, during the middle to late Miocene, the CAP had a paleoelevation between 1200 and 2000 m, a MAT of ~20 to 22C, and a MAP ranging from 550 to 1500 mm (fig. S5) (28, 29). Although other MAP estimates based on paleobotanical data from other CAP localities suggest lower precipitation ranges for the Miocene (30, 31), these values are never lower than modern records. The incongruence between the paleobotanical evidence from the middle to late Miocene and climate models for the CAP is a major question that needs to be addressed to predict the impact of future climate change in South America.

Our early Pliocene reconstructions indicate that, around 4.8 Ma, the NCAP had already attained near-modern elevations (~3800 m) and had a MAP slightly higher than modern (~850 mm). This higher-than-modern MAP is also supported by the fossil lacustrine biota (diatoms and ostracods) (15), stratigraphic (15), and isotopic data (8, 32). Given that ferns have higher water requirements than angiosperms (33), the high proportion of fern taxa present in Member C of the Descanso Basin compared with the modern puna (fig. S6B) offers additional evidence to support the high MAP estimations for the early Pliocene. The time spanned by Member C overlaps with the early Pliocene warm period (4.6 to 3.1 Ma). Our inferred wetter conditions for the NCAP during the early Pliocene contradict some global climate models that report permanent El Niolike conditions during the early Pliocene warm periodwhich implies lower-than-modern precipitation for the CAP(34, 35). However, whether these wetter conditions correspond to a large scale or a permanent state for the Altiplano during the Pliocene is unknown as numerous high-magnitude regional sea surface temperature (SST) excursions (36) and short-term climate perturbations (37) have also been reported for the Pliocene.

Separately, when comparing Pleistocene with Pliocene climate data, evidence from Pleistocene cores in the Altiplano suggest that glacial (interglacial) periods were wetter (drier) than modern (3840). Considering that SSTs for the Peru margin show that the early Pliocene was 3.3C warmer than the late Pleistocene (36) and that the elevation of the CAP has not changed substantially since the early Pliocene, paleoprecipitation would be expected to follow the interglacial pattern with drier regional conditions. However, our wetter CAP estimates for the early Pliocene suggest that this time period does not represent an extended interglacial period.

To address the similarities of paleofloras in members B and C to modern ecosystems, we used the Chaos dissimilarity index (CDI). We first calculate CDI for all extant sites in four elevation ranges (lowlands, lower montane, upper montane, and puna) to estimate the overall heterogeneity of the ecosystem, and then we compared the extant sites with the fossil palynological record from each member (Table 1). For Member B, the median CDI between palynoflora and extant lower montane sites was 0.8. This value is considerably higher than the median CDI of 0.3 obtained by all pairs of extant lower montane modern sites, suggesting high dissimilarity (Fig. 4A and Table 1). The median CDI between the Member B palynoflora and all extant upper montane sites was 0.7. This value differs substantially from the median CDI of 0.17 obtained by all pairs of extant upper montane modern sites, also suggesting high dissimilarity (Fig. 4B and Table 1). Thus, our results suggest that the palynoflora of Member B is not similar to extant Andean montane sites. This Miocene palynoflora is composed of common montane forest indicators such as Podocarpus, Hedyosmum, Bocconia, Mutisia, and Cyatheaceae, as well as higher elevation taxa such as Calamagrostis, Polylepis, and Valeriana. However, the oddity of this flora is emphasized by the co-occurrence of plants with traits of lowland ecosystems (e.g., 30-m legume trees and palms). Therefore, we propose that the Member B paleoflora does not have a modern analog and that it was more heterogeneous than modern montane ecosystems.

Floristic comparison of pollen records from members B and C and modern pollen records from sites located along the Central and Northern Andes (Fig. 1A). The altitudinal gradient includes sites from lowland, lower montane, upper montane, and puna ecosystems. The comparison was done using the CDI (0 = least dissimilar; 1 = most dissimilar).

Comparison of palynofloras from members B and C with extant pollen data from South America. The density plots show the probability distribution of the CDI (0 = least dissimilar; 1 = most dissimilar). Green curves represent CDI for within-biome comparison for modern altitudinal ranges (lower montane, upper montane, and puna). (A) Comparison for lower montane sites: Gray curve represents CDI between the fossil assemblage from the Miocene Member B and the modern lower montane sites. (B) Comparison for upper montane sites: Gray curve represents CDI between the fossil assemblage from the Miocene Member B and the modern upper montane sites. (C) Comparison for puna sites: Red curve represents CDI between the fossil assemblage from the Pliocene Member C and the modern puna sites.

For Member C, the median CDI between the palynoflora and all extant sites from the Puna was 0.2. This value was similar to the mean CDI of 0.16 obtained by all pairs of puna modern sites, suggesting low dissimilarity (Fig. 4C and Table 1). The macrofossil record of Member C also indicates a puna biome by the combined presence of Polylepis, Berberis, Ribes, and Polystichum. The Member C flora, therefore, represents the earliest record of the Puna, appearing shortly after the landscape is inferred to have reached near-modern elevations. This early puna differs from the modern puna by having a higher diversity and abundance of ferns (fig. S6B).

Late Miocene uplift shaped modern ecology: The rise of the Tibetan Plateau gave rise to the Asiatic monsoon system (41), while the uplift of the Andes created the South American Summer Monsoon (6). Moisture-laden winds flowing across Amazonia started to be deflected southward by the rising Andes between the Miocene and the Pliocene, drying out the highlands. This Andean uplift created new Neotropical niches for highland and dryland (Peruvian coastal rain shadow) species, while markedly reducing the area available to mid-montane species. Our data indicate that there would also have been a substantial reduction in potential biomass as forests supporting large trees were replaced by grasslands. The fact that current climate models underestimate Miocene precipitation so severely leads to underestimates of carbon storage and overestimates of the grasslands available to Miocene megafauna. Capturing the topographic effects of the Andes on atmospheric circulation accurately is essential to modeling both past and future climate change.

The coexistence approach proposed here offers a novel way to analyze macrofossil and palynological data and improves reconstructions of past climates and ecosystems. Our paleoelevation estimations calculated on the basis of the Neogene plant fossil record of the NCAP support evidence from previous studies based on isotopic data (8, 9) that proposed rapid surface uplift of approximately 2500 m in elevation between ~9.1 and 4.8 Ma. In contrast, our paleoprecipitation estimations strongly differ from those predicted by regional and global precipitation simulations for the Miocene and Pliocene [e.g., (3, 5, 29, 38)], revealing a marked precipitation decrease over the NCAP with the Andean uplift. Understanding the origin of these discrepancies can help reveal the driving forces controlling the climate of the Altiplano and ultimately the overall climate of South America. The uplift of the CAP shaped the biome distribution and climate of the region: Our reconstructions of a montane forest ecosystem without modern analogs in the Miocene and of a puna-like ecosystem in the early Pliocene highlight the dynamic nature of high elevation ecosystems and how these can change and evolve as a response to extreme modifications in elevation and climate.

Fieldwork. The Descanso-Yauri Basin, southern Peru, was explored during three field trips done during the dry season of years 2014, 2015, and 2016. We collected as many samples as possible during those field trips. All the samples were collected from members B and C. During the collection process, we kept a stratigraphic control by visiting the localities previously described by Kar et al. (8) and then by correlating stratigraphically the new localities found. We collected palynological samples, tuffs, and compressions and impressions of fossil leaves and fertile parts and permineralized wood. To correlate stratigraphically two of the new localities, we used radiometric dating of volcanic horizons (tuffs) present near the fossiliferous region.

Zircon U-Pb geochronology methods. Two tuff samples with abundant pumaceous fragments were collected to correlate two wood localities of the Colpamayo and San Miguel regions using U-Pb geochronology (Fig. 2). The sample from the Colpamayo region (STRI-MUSM 44449) was collected from a 60-cm conglomerate rich in pumaceous fragments, and the sample from the San Miguel region (STRI-MUSM 44448) was collected from a 50-cm-thick ashfall tuff. Zircon separates were performed by standard magnetic and gravity-based methods at the University of Rochester using a Frantz LB-1 separator and methylene iodide. Before U-Pb analysis by laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS), zircons were mounted in 1 epoxy plugs and characterized by cathodoluminescence (CL) imaging using a JEOL JSM-7100 electron microscope equipped with a field-emission gun and a Deben Centaurus CL detector at the Mackay Microbeam Laboratory, University of Nevada, Reno. All U-Pb isotopic measurements were performed by LA-ICP-MS at the Arizona Laserchron Center, University of Arizona using a Photon Machines Analyte-G2 Excimer laser system coupled to a Thermo Finnigan Element2 single collector ICP-MS. Analytical methods for the U-Pb analyses are outlined in Ibez-Mejia et al. (42) and Pullen et al. (43) (see the Supplementary Materials for detailed methods and reference material).

Palynological data. Palynological sample preparation included digestion of 10 g of rock in hydrochloric acid for 12 hours, addition of water, and decantation after 12 hours, followed by digestion in hydrofluoric acid for at least 24 hours, addition of water, and decantation of the acid solution after 24 hours. Sieving of the dissolved mineral fraction was initially done with a 250-m mesh to eliminate the thick fraction, followed by sieving through to 10-m mesh. Panning of the >10-m fraction was done in an ultrasonic equipment for recovery of the less dense organic matter fraction. This residue was cleaned in the ultrasonic equipment for some seconds, and the organic residue was concentrated by centrifugation, followed by mounting of a first cover slide in a solution of polyvinyl alcohol. Mounting of a second cover slide was done after oxidation with nitric acid, sealing with Canada balsam. Samples were processed at the laboratory of Paleoflora, Bucaramanga, Colombia. Palynomorphs were counted aiming for a minimum of 300 specimens per slide. Transects of all slides were made at 40 magnification using an Olympus BH-2 binocular scope to identify all pollen types. Bright-field microphotographs were obtained at 100 magnification using a Pixera Camera System attached to the Olympus scope.

All the taxonomic affinities proposed for the palynomorphs are reported in alphabetical order by family and then by genus (file S5). The proposed botanical affinities at the family, genus, and species level were proposed by comparing the fossil taxa with pollen atlases of extant taxa of the region and the pollen reference collection of A. Graham at the Center for Tropical Paleoecology and Archaeology, from the Smithsonian Tropical Research Institute, Panama (see the Supplementary Materials for reference material). The taxonomic status of botanical names was updated consulting the Tropicos database. Some associated organisms and those undetermined palynomorphs are included as assemblage of possible indicators of environmental conditions.

Macrofossil data. Macrofossil compressions and impressions were studied at the Paleontological Collection of the Museo de Historia Natural de la Universidad Nacional Mayor de San Marcos, Lima, Peru (DPV-MHN-UNMSM). Sediment was carefully removed using an air scriber to expose possible attachments and the maximum number of features. Specimens were observed using a Leica EZ4 HD coupled to an integrated camera of 5.0-megapixel with a complementary metal-oxide semiconductor sensor. Photographs were taken with varied low-angle light. Each leaf morphotype was described following the terminology of the Manual of Leaf Architecture (44). Comparisons with extant taxa were made with herbarium material from the San Marcos Herbarium (MHN-UNMSM, Lima, Peru), and virtual collections were accessed through the JSTOR Global Plants database and literature.

Permineralized wood was initially cut with diamond saws at the Museum of Natural History of Universidad de San Marcos of Lima (MUSM). Thin sections were prepared following the standard techniques described in detail by Boonchai (45). The anatomy of the samples was described following the International Association of Wood Anatomists (IAWA) list of microscopic features for hardwood identification (46). Preliminary comparisons with modern and fossil wood were done using the InsideWood database (47), and additional comparisons were done on the basis of literature.

Using linear regression models that account for the allometric relation between stem diameter and tree height in modern trees, we estimated the approximate height of one complete fossil sample (48, 49). Four different regression models were used for the calculation, one that uses a global dataset (48), while the other three used large tropical datasets from west Amazonia, tropical South America, and the pantropical region (49).

Macrofossil compressions and impressions and permineralized wood were organized and described by morphotypes following the method described by Peppe et al. (50). Each morphotype has a three-letter prefix (DSB or DSC) based on the name of the formation and the member from which they belong to, plus a number starting from one. A systematic affinity was proposed for each morphotype. Species names were not proposed for any of the morphotypes described because this required extensive research into the nomenclature of the taxon, previous fossil descriptions, and phylogenetic relationships (50), which were topics out of the scope of this paper.

A new coexistence method was developed here to estimate paleoelevation and paleoprecipitation using the macrofossil and palynological data collected in this study. Climate and elevation information associated with the nearest living relatives of the identified palynomorphs and macrofossil taxa was used to estimate a mutual climate and altitudinal range. Two distribution datasets (palynological and macrofossil) were described through bivariate probability density distributions of modern taxon occurrence both of our variables of interest.

Given that palynological counts were uneven because of differential preservation of samples and also aiming to make palynology and macrofossil estimates comparable, palynological relative abundance was not considered representative of environmental conditions, and the dataset was therefore transformed to presence-absence. Palynological samples were filtered to consider only samples that met the following criteria: (i) The sample was composed of taxa that were represented in the modern pollen and spore dataset, (ii) at least three taxa were represented in the sample, and (iii) at least one of the taxa present in the sample was identified at the genus level. For macrofossil material, taxonomic identifications at the family and genus level were used for the paleoclimatic analysis, and in instances of taxonomic uncertainty in the identifications, the probability density distributions were based on the modern distribution of multiple taxa that were phylogenetically related.

The distribution of taxa found in the fossil datasets (palynology and macrofossils) was described through bivariate probability density distributions (PDD) of modern taxon occurrence along elevation and precipitation gradients. These distributions were described through a bivariate Gaussian kernel density estimator (51) defined byf(x;H)=1ni=1nKH(xXi)where f(x;H) defines the probability density at a point x = (x1, x2)T; Xi = (Xi1, Xi2)T is the vector of realizations at points i = 1,2, , n, for variables X1 (elevation) and X2 (precipitation).

H=|h100h2| is a matrix of bandwidths; in our case, h1 = h2 = 250.

K(x)=12e12xTx is the Gaussian kernel smoother.

Georeferenced occurrences of the nearest living relatives of both palynological and macrofossil records were extracted from the BIEN R package (52) and were described in terms of the elevation and MAP at the occurrence points (data from WordClim) (53). Thus, two reference libraries of bivariate environmental distributions were built: one for palynology (PDFpaly) and one for macrofossils (PDFmacro). In instances of taxonomic uncertainty in the identification of macrofossils, PDFmacro was based on the modern distribution of multiple taxa that were phylogenetically related. The paleoenvironmental reconstruction for a given fossil palynology sample consisted of PDFfossil palyj derived from a nonweighted finite mixture (54) of the modern PDFpaly of sj pollen and spore taxa present in the fossil sample j. Individual sample estimates from the same stratigraphic member (T) were, in turn, mixed to obtain a PDFTPDFfossil palyj=k=1sjPDFpalykleadsPDFT=j=1nTPDFfossil palyjwhere PDFfossil palyj is the probability density function (PDF) of sample j composed of sj pollen taxa; PDFpalyk is the PDF of modern pollen taxon k with k = 1,2, , s.; PDFT is the PDF of stratigraphic member T, represented by nT samples.

The paleoenvironmental reconstruction based on macrofossils consisted of a single PDF for each stratigraphic member (PDFmacroB and PDFmacroC), as followsPDFmacro=i=1sPDFmacrokwhere PDFmacrok is the PDF of modern pollen taxon k with k = 1,2, , s. PDFpaly and PDFmacro of each sample represent the probability distribution of data pairs of precipitation and elevation conditioned to the presence of a given set of taxa s. The credible environmental area within each PDF was obtained by trimming the margins of the bivariate surface, keeping the area of the distribution associated with a cumulative probability higher than 0.5. On the basis of the probability densities of this area, the marginal distribution of elevation and precipitation for each sample was extracted, obtaining one macrofossil-based marginal distribution per variable per stratigraphic member and nB and nC palynology-based marginal distributions per variable. For summarizing palynology-based estimates, marginal distributions were subsequently mixed.

Marginal distributions of elevation estimates based on palynology and macrofossils were corrected taking into account the global temperature within the estimated age range for each stratigraphic member between 12.87 and 8.4 Ma for Member B and 5.3 and 4.3 Ma for Member C, global temperature anomaly data from Hansen et al. (55). Observed global temperature anomalies within the time period of each stratigraphic member were randomly sampled, converted to elevation using a lapse rate of 0.6C (56), and added to our estimates (global temperature was warmer).

Final estimates of both elevation and precipitation per stratigraphic member resulted from mixing palynology- and macrofossil-based marginal distributions. Modern estimates were based on mixing palynology-based, cross-validated estimates with distributions of modern climatic data (this latter step was necessary because of the lack of macrofossils for the Holocene). Such mixture resulted from an addition of individual densities weighted by a proportion factor, as followsMDparameter=1*MDpaly+2*MDmacrowhere MDparameter is the marginal distribution of temperature or precipitation; MDpaly is the pollen-based marginal distribution of temperature or precipitation; MDmacro is the macrofossil-based marginal distribution of temperature or precipitation; 1 and 2 are the representation factors for palynology and macrofossil data, respectively, with 1 + 2 = 1.

This mixture model was developed to evaluate the proportion or weight that each proxy (palynological versus macrofossil) should be given to better estimate paleoclimatic conditions because this information was unknown from the literature. 1 was set to vary from zero (estimation based exclusively on macrofossils) to one (estimation based exclusively on palynology) aiming to evaluate the stability of estimates. When estimates varied monotonically with , we selected a solution based on a completely proportional mixture (i.e., 2= 1= 0.5), whereas when the estimation was notoriously affected by the selection of , the inflection point was preferred. The environmental space for each time period was represented as the conjunction of the interquartile range of elevation and precipitation.

Fossil palynological data from members B and C were used for a comparison with modern pollen data. The modern pollen dataset consisted of palynological counts grouped by family from South American sites that span an elevation gradient from 100 to 4700 m. We used the CDI (57) to estimate similarities between modern data and modern versus fossil palynological data. The parameters used to estimate the CDI here were the asymmetric Sorensen-type using probability estimations and taxa abundance data. This analysis was implemented in R using the CommEcol package (58).

Because the modern pollen dataset did not include spore counts, we also did a comparison of the proportion of angiosperms and ferns present in the palynological record from Member C and modern data from the Global Biodiversity Information Facility (59) using relative abundances. Modern data were extracted from a polygon that included the Altiplano region and the eastern and western flanks of the Andes. The polygon was between 17 and 11 of latitude and 74 and 69 of longitude. Within the analysis, a quaternary sample that was collected in the Espinar region was also included in the comparison. Last, the inferred habit from the taxa represented in the fossil palynological samples was also compared to evaluate changes in the ecosystems.

Acknowledgments: We thank R. Salas-Gismondi, J. Tejada, and M. Arakaki for providing help in the field and space at the Museum of Natural History of Lima, Peru; G. Umasi for valuable help in the field and the Espinar community for welcoming us; D. Caballero for help with a palynological analysis; N. Jud for helpful discussions about the wood samples; T. Jordan, M. A. Gandolfo, and K. Nixon for providing guidance and reviewing the manuscript multiple times; N. Kar for providing geological information for the region; J. Houston for providing photo rights for the cover photo; and C. Hoorn and an anonymous reviewer for insightful comments on this manuscript. Funding: C.M. thanks the Fulbright-Colciencias Doctoral Fellowship, from Colombia; the College of Agriculture and Life Sciences, the Plant Biology Section, the Harold E. Moore, Jr. Memorial and Endowment Funds, and the Cornell Graduate School Research Travel Grant, Cornell University, USA, for providing facilities, travel and research support. C.J. was funded by the Smithsonian Institution, the National Geographic Society (NGS 2014 9537-14), the Anders Foundation, the 1923 Fund, and the Gregory D. and Jennifer Walston Johnson. A.C.-M. was supported by Programa de Apoyos para la Superacin del Personal Acadmico de la UNAM (PASPA, DGAPA). M.I.-M. acknowledges funding by NSF-EAR 1926124. Author contributions: C.M. and C.J. led the writing and had significant contributions from W.C., M.B.B., and M.I.-M. C.M., and C.J., designed the study and organized and conducted fieldwork. W.C. contributed to project design. A.C.-M. and C.M. contributed to the paleoprecipitation and paleoelevation analyses. C.M. described and analyzed wood fossil samples. J.E.M. and C.J. contributed to palynological analyses. F.M., C.J., and C.M. contributed to the stratigraphic analyses and fieldwork. A.A. and C.M. prepared and described the leaf macrofossils. M.I.-M. and F.M. conducted the zircon U-Pb geochronologic analyses. M.B., A.C.-M., C.M., and C.J. contributed to the comparison with modern pollen. C.J and C.M. provided financial support. All authors contributed to the manuscript and figure preparation. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested to the authors.

Excerpt from:
Neogene precipitation, vegetation, and elevation history of the Central Andean Plateau - Science Advances