# REVISITING THE BIOME CONCEPT WITH A FUNCTIONAL LENS

EDITED BY : Daniel M. Griffith, Christopher J. Still and Colin P. Osborne PUBLISHED IN : Frontiers in Ecology and Evolution

#### Frontiers Copyright Statement

© Copyright 2007-2019 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.

The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.

Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.

Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.

As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.

All copyright, and all rights therein, are protected by national and international copyright laws.

The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use. ISSN 1664-8714 ISBN 978-2-88945-930-8 DOI 10.3389/978-2-88945-930-8

### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

## Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# REVISITING THE BIOME CONCEPT WITH A FUNCTIONAL LENS

Topic Editors:

Daniel M. Griffith, Oregon State University, United States Christopher J. Still, Oregon State University, United States Colin P. Osborne, University of Sheffield, United Kingdom

Image: The Grumeti River, Tanzania (Google Earth Image 2019 CNES / Airbus).

Early biogeographers such as Alexander von Humboldt recognized the broad-scale coupling of vegetation and climate. This observation shaped the modern biome concept which organizes ecosystems by assumed relationships to environmental controls. This approach has been criticized for missing key impacts on the distribution and functioning of biomes like historical contingency, biogeographic history, disturbance ecology, and evolution. Are biomes still a convenient framework for organizing our understanding of biodiversity? What factors determine the functional differences among and within biomes, and at what spatial, temporal, and phylogenetic scales are those drivers most important? How can we better represent the functional characteristics and dynamics of ecosystems? This Research Topic highlights the latest discussions and research on biomes, drawing from a wide range of approaches spanning from macroecology and phylogeography to remote sensing and modelling ecosystem responses to global change.

Citation: Griffith, D. M., Still, C. J., Osborne, C. P., eds. (2019). Revisiting the Biome Concept with a Functional Lens. Lausanne: Frontiers Media. doi: 10.3389/978-2-88945-930-8

# Table of Contents

*06 Editorial: Revisiting the Biome Concept With a Functional Lens* Daniel M. Griffith, Christopher J. Still and Colin P. Osborne

## CHAPTER 1

## BIOME DISTRIBUTIONS AND FUNCTIONAL DIVERSITY

*09 Plant Functional Diversity and the Biogeography of Biomes in North and South America* Susy Echeverría-Londoño, Brian J. Enquist, Danilo M. Neves, Cyrille Violle,

Brad Boyle, Nathan J. B. Kraft, Brian S. Maitner, Brian McGill, Robert K. Peet, Brody Sandel, Stephen A. Smith, Jens-Christian Svenning, Susan K. Wiser and Andrew J. Kerkhoff


David L. Fox, Stephanie Pau, Lyla Taylor, Caroline A. E. Strömberg, Colin P. Osborne, Catherine Bradshaw, Stephen Conn, David J. Beerling and Christopher J. Still

## CHAPTER 2

## FUNCTIONAL VARIATION WITHIN AND AMONG BIOMES

*63 Grass Functional Traits Differentiate Forest and Savanna in the Madagascar Central Highlands*

Cédrique L. Solofondranohatra, Maria S. Vorontsova, Jan Hackel, Guillaume Besnard, Stuart Cable, Jenny Williams, Vololoniaina Jeannoda and Caroline E. R. Lehmann

*77 Functional Traits of Trees From Dry Deciduous "Forests" of Southern India Suggest Seasonal Drought and Fire are Important Drivers*

Jayashree Ratnam, S. K. Chengappa, Siddarth J. Machado, Nandita Nataraj, Anand M. Osuri and Mahesh Sankaran

*83 Inserting Tropical Dry Forests Into the Discussion on Biome Transitions in the Tropics*

Kyle G. Dexter, R. Toby Pennington, Ary T. Oliveira-Filho, Marcelo L. Bueno, Pedro L. Silva de Miranda and Danilo M. Neves

*90 Transplant Experiments Point to Fire Regime as Limiting Savanna Tree Distribution*

Nicola Stevens, Sally A. Archibald and William J. Bond

*103 Intraspecific Trait Variability in* Andropogon gerardii*, a Dominant Grass Species in the US Great Plains*

Seton Bachle, Daniel M. Griffith and Jesse B. Nippert

# Editorial: Revisiting the Biome Concept With A Functional Lens

Daniel M. Griffith<sup>1</sup> \*, Christopher J. Still <sup>1</sup> and Colin P. Osborne<sup>2</sup>

<sup>1</sup> Forest Ecosystems and Society, Oregon State University, Corvallis, OR, United States, <sup>2</sup> Department of Animal and Plant Sciences, University of Sheffield, Sheffield, United Kingdom

Keywords: biome, ecosystem function, functional traits, vegetation, evolutionary history

**Editorial on the Research Topic**

### **Revisiting the Biome Concept With A Functional Lens**

Early biogeographers such as Alexander von Humboldt recognized the broad-scale coupling of vegetation and climate (von Humboldt, 1806). This observation shaped the modern biome concept which organizes ecosystems by assumed relationships to environmental controls. Biomes are essential constructs for understanding vegetation distributions, the evolutionary patterns that shape species pools (Crisp et al., 2009; Cornwell et al., 2014), and the environmental impacts of human activities (Olson et al., 2001; Mucina, 2019), among other applications. However, ecologists recognize that there are regions, especially in the tropics, where vegetation may not deterministically relate to climate (Whittaker, 1975; Staver et al., 2011; Moncrieff et al., 2015). The biome concept is operationalized in practice as a static classification of the land surface. Process models rely on these classifications to summarize vegetation into Plant Functional Types (PFTs) which form the basis for representing ecosystem function and biogeochemical rates. Recently this approach has been criticized for missing key impacts on the distribution and functioning of biomes like historical contingency, biogeographic history, disturbance ecology, and evolution (reviewed in Higgins et al., 2016; Pausas and Bond, 2018). Thus, further research is required to better define biomes based on species composition and phylo-functional diversity, as well as better understand the drivers of biome boundaries and functioning within and among biomes.

A new understanding of biomes is crucial for appropriate prediction of future environmental change and global biogeochemical cycle modeling based on highly abstracted PFTs (Higgins et al., 2016; Still et al., 2018). In this issue, we present synthetic research ranging from continentalscale biogeography of biomes (e.g., Echeverría-Londoño et al.; Pinto-Ledezma et al.) to functional assessments of individual dominant species (e.g., Bachle et al.). These studies combine functional data with species distributions and phylogenies to provide new insight into the nature of biomes and how we can best capture the functional impacts of unique biogeographic histories (e.g., rare long-distance dispersal Deng et al.). They also indicate a major need for field data to fill gaps in datasets and for model parameterization. Focusing on all of North and South America, Echeverría-Londoño et al. compare and contrast functional diversity across biomes. They analyze the distributions of over 80,000 plant species combined with functional trait data for an ∼8,000 species subset. They report a general relationship between species range size and functional distinctiveness. Rare species are functionally distinct, whereas common species are functionally similar within each biome.

#### Edited and reviewed by:

Peter Convey, British Antarctic Survey (BAS), United Kingdom

> \*Correspondence: Daniel M. Griffith griffith.dan@gmail.com

#### Specialty section:

This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution

> Received: 27 March 2019 Accepted: 12 April 2019 Published: 01 May 2019

#### Citation:

Griffith DM, Still CJ and Osborne CP (2019) Editorial: Revisiting the Biome Concept With A Functional Lens. Front. Ecol. Evol. 7:144. doi: 10.3389/fevo.2019.00144

Extant ecosystem function is a product of the assembly processes that shape the structure of communities (Pennington et al., 2004; Higgins, 2017; Mucina, 2019) and the evolutionary processes that interact with filters and species relations (HilleRisLambers et al., 2012; Cavender-Bares et al., 2016). Here, Pinto-Ledezma et al. use a phylo-functional approach to partition beta diversity into two major sources of community compositional change: nestedness, representing change attributed to species loss; and turnover, arising from species replacement. The analysis spans the biomes of North America, and reveals that diversity in species-rich biomes with stable environments tends to arise from a combination of speciation processes and local environmental sorting. Together these produce species turnover along environmental gradients. Differences in biogeographic histories predict patterns of functional similarity among and within these biomes, with nestedness being more important for functional change than turnover.

Biome history is viewed through a different lens by Fox et al. who examine the geologic history of tropical grassy biomes. These biomes were assembled via the increasing abundance of grasses using the C<sup>4</sup> photosynthetic pathway, which came to dominate open, tropical habitats during the late Miocene (Osborne and Beerling, 2006; Edwards et al., 2010; Strömberg, 2011). This biome assembly caused major shifts to presentday tropical carbon cycling (Still et al., 2003) as well as past and present fire regimes (Scheiter et al., 2012). Yet the climate responses of disjunct savanna ecosystems differs considerably across continents (Lehmann et al., 2014). Fox et al. analyze more than 2,600 fossil isotope values to document global Neogene variation in C<sup>4</sup> grass abundance. They find significant but fairly weak agreement between isotope proxy values and climate-driven model predictions of varying complexity. This suggests historic roles for disturbance, biogeographic history, and local ecology in influencing patterns of ecosystem change and function, matching the situation in modern grassy ecosystems (Griffith et al., 2015). These factors are neither fully incorporated into current ecosystem models nor well represented in biome classifications.

The way in which biomes are classified has major impacts for conservation, landscape management, and projected ecosystem change (Pennington et al., 2004; Banda et al., 2016; Lehmann and Parr, 2016; Moncrieff et al., 2016; Griffith et al., 2017). Savannas are often misclassified as degraded forests, resulting in mismanagement of fire regimes and tree planting in ancient grassy ecosystems. Tropical and subtropical savannas are distinguishable from forests by a flammable C<sup>4</sup> grass understory and trees adapted to fire, despite overlap in tree cover values (Scholes and Archer, 1997; Ratnam et al., 2011). Here, Solofondranohatra et al. use grass phylogeny, vegetation surveys, and trait data to extend this concept to differentiate the evolutionary history and function of the understory. They provide new evidence to show that woodland regions of Madagascar are phylogenetically and functionally savannas and not degraded forests, as is asserted by many current biome classifications. In southern India, Ratnam et al. examined a range of vegetation types and found that much of the vegetation previously classified as forests was functionally savannas, having tree species with traits associated with frequent fires. This suggests that current fire suppression practices may be inappropriate for large areas of this region. Conversely, Dexter et al. show that dry forests in South America are functionally distinct from both savannas and moist forests. Dry forests are characterized by tree species adapted to seasonal drought stress and high soil fertility, but lack fire as a significant ecological driver. This means that conservation of dry forest ecosystems also requires unique management practices. Furthermore, these functional differences have a major influence on the nature of biome boundaries and transitions among biome states. For example, Dexter et al. suggest that transition zones between dry and moist forests may be dominated by changes in water availability, whereas transitions between savannas and moist forests may be sharp boundaries characterized by feedbacks and alternative stable states (Hoffmann et al., 2012).

Phylo-functional variation within biomes likely influences the response of ecosystems to climate change and helps explain different functioning among the same biome type on different continents (Lehmann et al., 2014). A major source of variation in savannas across continents is the species pool of trees from which fire-tolerant savanna trees evolved. Stevens et al. provide evidence from a continental-scale transplant experiment that the range of two dominant African savanna trees is most limited by their ability to escape the fire trap. This work further highlights that models of species distributions in savannas, and models of savanna ecosystems in general, will be inappropriate when only considering climatic descriptors of species ranges. Intraspecific trait variation is another source of functional differentiation within biomes. Bachle et al. assembled available data for key functional traits of Andropogon gerardii, a dominant grass species of the US Great Plains. Their synthesis suggests that this grass's high abundance and widespread distribution is enabled by its functional attributes and the potential for intraspecific variation to buffer populations against climatic variation. This finding echoes the general finding of Echeverría-Londoño et al. that widespread species will be less specialized functionally. Furthermore, Bachle et al. expose a surprising dearth of functional data for A. gerardii, especially given that they focus one of the most well-studied grasses and ecosystems in the world. This reinforces the observation from data syntheses (e.g., Echeverría-Londoño et al.; Pinto-Ledezma et al.), and in fact all studies in this special issue, that vastly more phylogenetic and functional data are required to appropriately understand and project biome function and distribution into the future.

This special issue highlights that, across multiple levels of biological organization, biomes still provide an important conceptual framing of global ecology. Biome classifications are most successful when they include functional and phylogenetic information, and when they allow biomes and boundaries to emerge from species-level data and ecological interactions, rather than imposed by PFT associations with climate and soils. Functional biome classification and mapping approaches will better encapsulate the evolutionary history that produced modern ecosystem function. These studies outline a massive challenge for ecologists and modelers in order to help predict and mitigate rapid and potentially irreversible modifications to the functioning of Earth's biomes.

## AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

## REFERENCES


## FUNDING

This work was supported by NSF award 1342703 (CS and DG).

## ACKNOWLEDGMENTS

Thank you to the contributing authors for making this special issue successful. Thank you to the Frontiers in Ecology and Evolution editorial staff and reviewers for their assistance.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Griffith, Still and Osborne. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Plant Functional Diversity and the Biogeography of Biomes in North and South America

Susy Echeverría-Londoño<sup>1</sup> \*, Brian J. Enquist <sup>2</sup> , Danilo M. Neves 2,3, Cyrille Violle<sup>4</sup> , Brad Boyle<sup>2</sup> , Nathan J. B. Kraft <sup>5</sup> , Brian S. Maitner <sup>2</sup> , Brian McGill <sup>6</sup> , Robert K. Peet <sup>7</sup> , Brody Sandel <sup>8</sup> , Stephen A. Smith<sup>9</sup> , Jens-Christian Svenning10,11, Susan K. Wiser <sup>12</sup> and Andrew J. Kerkhoff <sup>1</sup>

<sup>1</sup> Department of Biology, Kenyon College, Gambier, OH, United States, <sup>2</sup> Department of Ecology and Evolutionary Biology, University of Arizona, Tucson, AZ, United States, <sup>3</sup> Department of Botany, Federal University of Minas Gerais, Belo Horizonte, Brazil, <sup>4</sup> Centre d'Ecologie Fonctionnelle et Evolutive (UMR 5175), CNRS, Université de Montpellier, Université Paul Valéry, Montpellier, France, <sup>5</sup> Department of Ecology and Evolutionary Biology, University of California, Los Angeles, Los Angeles, CA, United States, <sup>6</sup> School of Biology and Ecology, University of Maine, Orono, ME, United States, <sup>7</sup> Department of Biology, University of North Carolina, Chapel Hill, NC, United States, <sup>8</sup> Department of Biology, Santa Clara University, Santa Clara, CA, United States, <sup>9</sup> Department of Ecology and Evolutionary Biology, University of Michigan, Ann Arbor, MI, United States, <sup>10</sup> Department of Bioscience, Center for Biodiversity Dynamics in a Changing World (BIOCHANGE), Aarhus University, Aarhus, Denmark, <sup>11</sup> Section for Ecoinformatics and Biodiversity, Department of Bioscience, Aarhus University, Aarhus, Denmark, <sup>12</sup> Landcare Research New Zealand, Lincoln, New Zealand

#### Edited by:

Daniel M. Griffith, Oregon State University, United States

#### Reviewed by:

Joaquín Hortal, Museo Nacional de Ciencias Naturales (MNCN), Spain Mark E. Olson, National Autonomous University of Mexico, Mexico

\*Correspondence:

Susy Echeverría-Londoño susyelo@gmail.com

#### Specialty section:

This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution

Received: 02 June 2018 Accepted: 03 December 2018 Published: 18 December 2018

#### Citation:

Echeverría-Londoño S, Enquist BJ, Neves DM, Violle C, Boyle B, Kraft NJB, Maitner BS, McGill B, Peet RK, Sandel B, Smith SA, Svenning J-C, Wiser SK and Kerkhoff AJ (2018) Plant Functional Diversity and the Biogeography of Biomes in North and South America. Front. Ecol. Evol. 6:219. doi: 10.3389/fevo.2018.00219 The concept of the biome has a long history dating back to Carl Ludwig Willdenow and Alexander von Humboldt. However, while the association between climate and the structure and diversity of vegetation has a long history, scientists have only recently begun to develop a more synthetic understanding of biomes based on the evolution of plant diversity, function, and community assembly. At the broadest scales, climate filters species based on their functional attributes, and the resulting functional differences in dominant vegetation among biomes are important to modeling the global carbon cycle and the functioning of the Earth system. Nevertheless, across biomes, plant species have been shown to occupy a common set of global functional "spectra", reflecting variation in overall plant size, leaf economics, and hydraulics. Still, comprehensive measures of functional diversity and assessments of functional similarity have not been compared across biomes at continental to global scales. Here, we examine distributions of functional diversity of plant species across the biomes of North and South America, based on distributional information for > 80,000 vascular plant species and functional trait data for ca. 8,000 of those species. First, we show that despite progress in data integration and synthesis, significant knowledge shortfalls persist that limit our ability to quantify the functional biodiversity of biomes. Second, our analyses of the available data show that all the biomes in North and South America share a common pattern–most geographically common, widespread species in any biome tend to be functionally similar whereas the most functionally distinctive species are restricted in their distribution. Third, when only the widespread and functionally similar species in each biome are considered, biomes can be more readily distinguished functionally, and patterns of dissimilarity between biomes appear to reflect a correspondence between climate and

**9**

functional niche space. Taken together, our results suggest that while the study of the functional diversity of biomes is still in its formative stages, further development of the field will yield insights linking evolution, biogeography, community assembly, and ecosystem function.

Keywords: biogeography, biomes, functional traits, hypervolumes, macroecology, plant functional diversity

## INTRODUCTION

Ecologists and biogeographers organize land plant biodiversity into climatically-determined biomes, with physiognomies characterized by the growth forms and functional traits of the dominant species (Moncrieff et al., 2016). Indeed, the biome concept has a long history dating back to Carl Ludwig Willdenow and Alexander von Humboldt. Willdenow recognized that similar climates support similar vegetation forms, and Humboldt observed the widespread association between plant distribution, physiognomy, and environmental factors. Overall, the biome concept reflects the assumption that similar environmental pressures select for species with similar functional attributes, independent of their evolutionary history. At the same time, Earth's extant biomes are to some extent phylogenetically distinct, with many/most of the characteristic species being drawn from specific lineages exhibiting key adaptations not just to climatic selection but to additional pressures like fire or megaherbivores (Woodward et al., 2004; Pennington et al., 2006; Donoghue and Edwards, 2014). Because biomes represent broad-scale regularities in Earth's vegetation, understanding functional differences among biomes is critically important to modeling the global carbon cycle and the functioning of the Earth system, including responses to anthropogenic global change (Bonan et al., 2012; van Bodegom et al., 2014; Xia et al., 2015).

The past two decades has seen rapid growth in understanding the functional dimension of plant biodiversity across continental scales (Swenson et al., 2012; Lamanna et al., 2014; Šímová et al., 2018). Extensive data synthesis and analysis efforts have defined a limited functional trait space incorporating variation in overall plant size, leaf economics, and hydraulics, known as functional trait spectra (Wright et al., 2004; Díaz et al., 2016). Functional trait spectra define physiological and ecological trade-offs that determine plant life-history strategies (Westoby, 1998; Reich et al., 2003; Craine, 2009) and their influence on community assembly (Kraft et al., 2008), ecosystem function (Lavorel and Garnier, 2002; Garnier et al., 2004; Kerkhoff et al., 2005; Cornwell et al., 2008), and even rates of evolution (Smith and Donoghue, 2008; Smith and Beaulieu, 2009). The generality and utility of plant functional trait spectra has hastened their incorporation into models of the global distribution of vegetation (van Bodegom et al., 2014), biogeochemical cycling (Bonan et al., 2012), and ecosystem services (Díaz et al., 2007; Cadotte et al., 2011; Lavorel et al., 2011). However, few studies have examined comprehensive measures of functional diversity at the biome scale across both dominant and subordinate species and life forms.

There are two contrasting sets of predictions about how functional diversity varies among climatically and physiognomically distinct biomes. On the one hand, if the biodiversity of biomes reflects variation in the available ecological niche space, less taxonomically diverse biomes should represent a smaller, largely nested subset of the functional space occupied by more diverse biomes. On the other hand, due to the global nature of fundamental trade-offs in plant structure and function, and similar selection pressures acting similarly across species, assemblages occupying different environments may in fact share similar areas of trait space (Reich et al., 2003; Wright et al., 2004; Díaz et al., 2016). Recent studies confirm that the relationships between climate and functional and taxonomic diversity (i.e., species richness) are complex and scale-dependent. In some cases, functional diversity closely tracks climatic gradients, with lower functional diversity in more variable and extreme conditions (Swenson et al., 2012; de la Riva et al., 2018), and the volume of functional space in local assemblages expands with increasing taxonomic richness (Lamanna et al., 2014; Li et al., 2018). However, analyses of regional-scale species pools suggest that large changes in species richness may be associated with minimal impacts on functional diversity (Lamanna et al., 2014; Šímová et al., 2015). Moreover, the responses of taxonomic and functional diversity to climate may differ among major plant growth forms (Šímová et al., 2018), and the relative diversity of different growth forms may change substantially along the climate gradients that define different biomes (Engemann et al., 2016).

Progress in analyzing functional diversity on continental to global scales has been limited in part by data shortfalls in cataloging the taxonomic, distributional, phylogenetic, and functional aspects of biodiversity (Hortal et al., 2015). Furthermore, even given the limited data available, significant informatics challenges are associated with standardizing and integrating large, disparate datasets describing the geographic distributions, functional traits, and phylogenetic relationships of species (Violle et al., 2014). Here we examine the distribution of functional diversity of plant species across the biomes of North and South America using the Botanical Information and Ecology Network database (BIEN; Enquist et al., 2016; Maitner et al., 2018), which assembles distributional and functional trait information on >100,000 species of land plants. Specifically, we examine variation in the functional diversity and distinctiveness of biomes, based on land plant species distribution maps at a scale of 100 × 100 km grid cells and a comprehensive dataset of six functional traits that reflect the major axes of variation in ecological strategies.

Our goals in this study are 3-fold. First, in order to highlight persistent data shortfalls, we document the extent of the available data characterizing the functional diversity and distinctiveness of biomes. Second, given the data available, we characterize the distribution of functional diversity within biomes across both dominant and subordinate growth forms. These analyses allow us to better quantify the functional distinctiveness of a biome by identifying the most common functional strategies of the most widespread species within it. Third, we ask whether biomes are in fact characterized by functionally distinct collections of species, based on the overlap of multidimensional hypervolumes in functional trait space.

## METHODS

## Distribution Data and Biome Classification

To reduce the effects of sampling bias characteristic of occurrence datasets compiled from multiple resources, we used the BIEN 2.0 range maps for 88,417 of plant species distributed in North and South America (Goldsmith et al., 2016). The BIEN database integrates standardized plant observations stemming from herbarium specimens and vegetation plot inventories. Species range maps in the BIEN databases were estimated using one of three approaches, depending on the occurrences available for each species. For species with only one or two occurrences (c. 35% of the species), the geographic range was defined as a square 75,000 km<sup>2</sup> area surrounding each data point. The geographic ranges for species with three or four occurrences were identified using a convex hull (c. 15% of the species). Finally, the range maps for species with at least five occurrences were obtained using the Maxent species distribution modeling algorithm, using 19 climatic layers as predictor variables and 19 spatial eigenvectors as filters to constrain over predictions by the models (see Goldsmith et al., 2016 for further details on range map methodology).

We overlaid the BIEN 2.0 plant species maps on a 100 × 100 km grid map with a Lambert Azimuthal Equal Area projection to obtain a presence/absence matrix of species for each grid cell. Based on the Olson et al. (2001) biome classification, we assigned each matrix grid cell to one of the biomes categories. Because of computational limitations for the subsequent analyses, we joined some biomes based on their similarities in climate and vegetation and literature to obtain a broad classification as described in **Table 1** (excluding Inland Water, Rock and Ice and Mangroves). The Chaco and Caatinga ecoregions were classified as Xeric Woodlands and Dry forest, respectively, following Prado and Gibbs (1993), Pennington et al. (2000), Banda et al. (2016), Silva de Miranda et al. (2018).

## Species Composition Among Biomes

To understand the variation in functional trait space among biomes, we first explored the differences in plant species composition among them. Based on the species list for each grid cell and biome, we defined each biome's characteristic species as those that have the maximum proportion of their range in that biome. These lists of geographically predominant species in a given biome were compared to the list of species in other biomes. This pairwise comparison provides a simple way to infer a directional taxonomic overlap among biomes (i.e., the TABLE 1 | Overview of the biome classification adopted in this study and the equivalent (Olson et al., 2001) biomes classification (excluding Inland Water, Rock and Ice, and Mangroves).


Because of computer limitations, some biomes were joined based on their climatic and vegetation similarities. The Chaco and Caatinga ecoregions were classified as Xeric Woodlands and Dry forest, respectively, following Prado and Gibbs (1993), Pennington et al. (2000), Banda et al. (2016), Silva de Miranda et al. (2018).

proportion of predominant species of a biome shared with another biome). Of the 88,417 species with available range maps, 44,899 species have ranges spanning more than one biome, and 43,518 species are endemic to a specific biome.

## Trait Data

We extracted all the trait information for plant species in the New World available in the BIEN 3.0 dataset (retrieved on 7 February 2018) giving a total of 80,405 species levels trait observations. We then filtered the information for six functional traits: maximum plant height (m), seed mass (mg), wood density (mg/cm<sup>3</sup> ), specific-leaf-area SLA (cm<sup>2</sup> /g) and leaf phosphorus and leaf nitrogen concentration per unit mass (mg/g). A total of 18,192 species-level observations were left after filtering. From these, the subset of 8,820 species with both range maps and trait information was used for further analyses.

To estimate the functional trait space for each biome, we required complete trait data. For this reason, we phylogenetically imputed missing trait data (Bruggeman et al., 2009; Penone et al., 2014; Swenson, 2014; Swenson et al., 2017) using the R package "Rphylopars" v 0.2.9 (Goolsby et al., 2017) and the recently published phylogeny of seed plants by Smith and Brown (2018) as a baseline. Phylogenetic imputation is a tool for predicting missing data in functional traits datasets based on the assumption that closely related species tend to have similar trait values (Swenson, 2014). We used the ALLBM tree (i.e., gene bank and Open Tree of life taxa with a backbone provided by Magallón et al., 2015) because it maximized the overlap between species with available trait and distributional information. After the trait imputation, a total of 7,842 species with complete trait information and range maps remained for further analysis.

## Trait Hypervolumes Measurement

We measured the relative functional diversity of biomes by calculating trait hypervolumes from species pools within biomes cells. Due to computational limitations, we constructed the trait hypervolumes using a random sample of 20% of the cells in each biome. The hypervolumes for each grid cell were estimated using the extracted six functional traits and the R package "hypervolume" (Blonder et al., 2014, 2018), using the Gaussian KDE method with the default Silverman bandwidth estimator. Seed mass, height and wood density were log-transformed, whereas SLA was squared-root transformed. All traits were scaled and centered before the analysis. Hypervolumes are reported in units of standard deviations (sd) to the power to the number of traits used (i.e., sd<sup>6</sup> ).

## Functional Distinctiveness and Widespreadness

Because biomes are characterized by their dominant vegetation, we also examined the geographic commonness and functional distinctiveness of species within biomes. Using species range maps and functional trait information, we estimated the functional distinctiveness and widespreadness of each species in each biome, following the conceptual framework of functional rarity by Violle et al. (2017). Using the whole set of (measured and imputed) traits, we first measured the Euclidean distance in standardized trait space between all pair of species. We then calculated the functional distinctiveness (Di) for each species in each biome as the average functional distance of a species to the N other species within the biome species pool. Di was scaled between 0 and 1 within each biome with lower values representing species that are functionally common (redundant) and higher values representing species that are functionally distinctive, compared to the other species in the biome. We also estimate the geographic widespreadness (Wi) of a species in a biome, measured as the number of grid cells occupied by the species within a biome over the total number of cells in that biome. A value of 1 indicates that the focal species is present in all grid cells covered by the biome. Both functional distinctiveness and geographic widespreadness were calculated using the R package "funrar" (Grenié et al., 2017).

Because we lack comprehensive measures of dominance based on local abundance or biomass, we used the measures Di and Wi to identify for each biome the most "common" species as those that are both geographically widespread (Wi > 0.5) and functionally similar (Di < 0.25). This last threshold was used since the third quantile of functional distinctiveness values among biomes ranged from 0.2 to 0.3. To discriminate functionally distinct species vs. functionally similar species in subsequent analyses, we, therefore, used a value of 0.25 as a cut-off (i.e., species with distinctiveness values < 0.25 were considered functional common or redundant within a specific biome).

## Functional Space and Biome Similarity

To estimate the overlap among biomes' hypervolumes we employed the Sørensen similarity index using (i) the total number of species and (ii) the list of species considered as functionally common and geographically widespread for each biome. Functional trait space similarity was calculated as the pairwise fractional overlap of hypervolumes among biomes. The fractional overlap was calculated by dividing twice the volume of the intersection of two hypervolumes by the volume of their unions. All hypervolumes were estimated using the R package "hypervolume" (Blonder et al., 2014, 2018), using the Gaussian KDE method with the default Silverman bandwidth estimator.

## RESULTS

We found substantial overlap in plant species composition across biomes, based on the intersection of modeled species ranges (see **Figure 1**, and **Supplementary Table 1**). This taxonomic overlap is greater within tropical and temperate biomes, with relatively few species being shared between these two climatic zones. Interestingly, xeric woodlands share species with both tropical and temperate biomes. Despite the high proportion of species characteristic of tropical moist forest (83%), this biome also shares large numbers of species with other tropical biomes such as dry forests, savannas, xeric woodlands and tropical grasslands. The high taxonomic overlap of temperate biomes such as temperate grasslands, Mediterranean woodlands, taiga and tundra suggests the potential for low functional distinctiveness, with some temperate biomes representing a poorer subset of the more species-rich, functionally-diverse tropical biomes.

With fewer than 10% of the mapped species represented, the trait data were quite sparse for our biome-scale species assemblages (**Figure 2A**). Moreover, the available trait data varied substantially across traits, plant clades, and biomes. Phylogenetically, leaf P and wood density were particularly poorly represented with whole clades lacking any available data. As a result, of the approximately 2/3 of the trait values for our 7,842 species that had to be imputed, many had to be assigned in the absence of any closely related taxa. Geographically, SLA, height, and seed mass are undersampled in the tropics and temperate South America, while leaf N and wood density are more evenly sampled among regions. Leaf P is poorly sampled across all biomes. Trait hypervolumes created from species pools of a random selection of 20% of cells for each biome mainly show differences between tropical and temperate/cold biomes (**Figure 2B**). Trait hypervolumes tend to be larger and more variable in tropical biomes, with the highest variation in moist and dry tropical forests. Among temperate/cold biomes, coniferous and temperate mixed forest support the largest trait hypervolumes. Interestingly, xeric shrublands exhibit a distribution of hypervolumes more similar to tropical than to temperate biomes.

The relationship between functional distinctiveness and geographic distribution is remarkably similar among biomes (**Figure 3**). In every biome, the vast majority of species are

both geographically restricted and functionally similar. At the same time, the most functionally distinctive species within every biome were generally geographically restricted as well, and the most geographically widespread species were almost always functionally similar.

We found considerable variation in the proportional distribution of growth forms within and among biomes. Woody species represent the dominant growth forms in tropical biomes in terms of species numbers, whereas herbaceous species dominate in temperate environments (**Figure 4**). When we considered only those species that are widespread and functionally common in each biome, the distribution of growth forms across biomes changed, especially in the proportion of trees, herbs and grasses. For instance, in tropical biomes, the proportion of trees decreased in each biome except for moist forest. In temperate biomes, the proportion of grasses increased, especially toward the polar regions. Interestingly, the distribution of growth forms in xeric woodlands more closely resembles the distribution in temperate biomes when only widespread and functionally common species are considered.

Pairwise comparisons of species composition among biomes reveal three main clusters representing the tropical, temperate and polar climatic zones (**Figure 5A**), reflecting the high taxonomic overlap within, but not between, these regions (**Figure 1**). Pairwise comparisons of trait hypervolumes among biomes show a less clear clustering of climatic zones (**Figure 5B**). In this case, the functional space of xeric woodlands overlaps significantly with temperate biomes. This analysis also reveals the overlap in functional space of taigas with temperate grasslands and mixed forests. The pairwise comparison of trait hypervolumes among biomes using only those species considered as functionally common and widespread shows less overlap in trait spaces among and within climatic zones (**Figure 5C**). However, even though xeric woodlands are now clustered with the rest of tropical biomes, these habitats along with tropical grassland exhibit great overlap in functional space with temperate biomes such as Mediterranean woodlands and temperate grasslands.

## DISCUSSION

Our analyses yield three important insights for understanding terrestrial biomes through a functional lens. First, we show that despite progress in the compilation and synthesis of primary biodiversity data, significant knowledge shortfalls persist that may limit our ability to quantify the functional biodiversity of biomes on continental to global scales. Second, our analyses of the available data nevertheless show that all of the biomes in North and South America share a remarkable common pattern in which the most geographically widespread species in any

proportion of species with trait information relative to the total number of species in the BIEN 2.0 database. (B) Distribution of trait hypervolumes of 20% of randomly selected 100 ×100 km cells in each biome. Hypervolumes are reported in units of standard deviations to the power of the number of traits used.

biome tend to share common functional traits while the most functionally distinctive species are invariably restricted in their distribution. Third, when only the widespread and functionally common species are considered, biomes can be more readily distinguished functionally, and patterns of dissimilarity between biomes appear to reflect a correspondence between climate and plant functional niche space. Taken together, our results suggest that while the study of the functional diversity of biomes is still in its formative stages, further development of the field will likely yield insights linking evolution, biogeography, community assembly, and ecosystem function.

## KNOWLEDGE SHORTFALLS

The BIEN database is specifically designed to close gaps in our knowledge of plant biodiversity, yet as Hortal et al. (2015) point out, interacting knowledge shortfalls lead to uncertainty in quantifying biodiversity at the largest scales. For example,

not only did fewer than 10% of our mapped species have both trait data and phylogenetic information available, but those missing trait data (∼66% of species' trait values) were quite structured, both phylogenetically and geographically. Thus, the Raunkiaerian shortfall in functional trait data (sensu Hortal et al., 2015) interacts with the Wallacean shortfall in distributional information and the Darwinian shortfall in phylogenetic understanding. And while we can use phylogenetic knowledge to reasonably impute those missing values (Swenson et al., 2017), several aspects of the resulting patterns that are important to biome classification remain uncertain. For example, the existing growth form data suggest that woody

species dominate in all tropical biomes, whereas the proportional diversity of herbaceous growth forms is much higher in temperate and polar biomes (see also Engemann et al., 2016; Zanne et al., 2018). Does this latitudinal increase in relative diversity of herbaceous plants reflect sampling biases and a lack of taxonomic knowledge about tropical herbaceous plants, or does it reflect the differing evolutionary and biogeographic histories in the Nearctic and Neotropical realms? Dominant growth forms are one of the key distinguishing features of biomes, so systematic changes in the most diverse growth forms (whether dominant or not) will necessarily influence our predictions about how functional diversity might change across biomes, and how those changes will affect ecosystem function and services.

We used species range maps as input data mainly to overcome the undeniable issue of sampling bias that is characteristic of datasets compiled from multiple sources (see **Supplementary Figure 1**). However, range maps drawn from species distribution models could represent an additional element of uncertainty in our results. For example, given the geographic resolution of this study and the spatial complexity of certain biomes, these models could have overestimated the geographical extent of some species from cells of a single biome to cells of nearby and interdigitated ones. Using occurrence data did not change the general results of this study (see **Supplementary Figures 2**–**6**); however, a slight decrease in the functional overlap among biomes could indeed reveal an overestimation of some species ranges (see **Supplementary Figures 2**, **5**, and **6**). Using occurrence information also changed the general relationships between functional distinctiveness and widespreadness in our results, due to the relative oversampling of temperate vs. tropical regions (see **Supplementary Figure 4**). In this case, our measure of widespreadness is limited by the extremely limited sampling of most species in tropical biomes. As datasets and sampling improve, there is likely to be scope to reduce these uncertainties in the future.

The other data gap that needs to be addressed to understand the functional diversity of biomes is the Prestonian shortfall in abundance information. Abundance is particularly important in studies of functional diversity because the traits that matter most for community assembly, ecosystem functioning, and biogeochemical cycling are those of the most dominant species. Because we based our assemblages on range maps derived from species distribution models, we could only address the geographic component of commonness, but the positive relationship between local abundance and geographic range, the so-called occupancy-abundance relationship, suggests that the most widespread species in each biome may frequently be among the more abundant as well (Borregaard and Rahbek, 2010; Novosolov et al., 2017). Because biomes are generally characterized by the dominant growth forms in the region, integrating abundance information might result in biomes showing less functional overlap. However, even though BIEN has compiled voluminous plot data that allow estimates of local abundance, making reasonable comparisons across different regions and growth forms is hampered by uneven sampling, regional differences in gamma diversity, and incommensurate methods of quantifying abundance.

## COMMON PATTERNS WITHIN BIOMES

Despite the persistent gaps in the available data, our analyses uncovered a previously undocumented relationship common to all the biomes of the Western Hemisphere: in any biome, the most widespread species also tend to exhibit low functional distinctiveness, whereas the most functionally distinctive species are almost invariably restricted in their distribution (**Figure 3**). This "occupancy-redundancy" relationship may suggest that the climatic conditions prevalent within a biome select for a set of common characteristics, with more functionally distinctive species being restricted to a rarer subset of habitats within the biome. The overall prevalence of functionally similar species in all biomes and the "occupancy-redundancy" relationship are both consistent with a recent global analysis of community level functional diversity that suggests that habitat filtering leads to the coexistence of functionally similar species (Li et al., 2018) as well as studies showing that functional redundancies increase community stability (e.g., Walker et al., 1999; Pillar et al., 2013). Moreover, together with the high degree of both taxonomic and functional overlap among biomes (**Figures 1**, **5**), the fact that common, widespread species are functionally similar reinforces the notion that land plants across a wide range of environmental conditions share common characteristics near the core of the functional trait spectrum (Wright et al., 2004; Díaz et al., 2016). Unlike Umaña et al. (2017), our results show that rare, geographically restricted species may or may not be functionally distinct from more widespread and common species. Thus, the question remains open whether more functionally distinctive species are specialists in particular environments or whether their trait combinations result in demographic or physiological trade-offs that limit their geographic distribution.

## COMPARISONS BETWEEN BIOMES

Our comparison of trait hypervolume distributions across biomes (**Figure 2B**) is consistent with the observation that more species-rich environments are also more functionally diverse (Swenson et al., 2012; Lamanna et al., 2014; Li et al., 2018; Šímová et al., 2018). Our results are more equivocal concerning the hypothesis that seasonal and extreme climatic environments consistently limit the functional diversity of species (de la Riva et al., 2018). All tropical biomes display higher average functional diversity than all temperate biomes, and the polar biomes do display the smallest hypervolumes. However, within each group, drier or more seasonally variable biomes do not always display smaller hypervolumes (e.g., dry forests, xeric woodlands), as we would have expected following the tolerance hypothesis by Currie et al. (2004), whereas optimal climatic conditions support more combinations of physiological parameters. From these results alone we cannot determine whether temperate and polar biomes are less taxonomically diverse because of limits on functional niche space, or whether their functional hypervolumes are small because they are not taxonomically diverse. Null-modeling approaches could potentially help to disentangle taxonomic and functional diversity (Swenson et al., 2012; Lamanna et al., 2014; Šímová et al., 2018), but such an

FIGURE 5 | (A) Pairwise dissimilarity in species composition among biomes. (B) Pairwise dissimilarity in trait hypervolumes (1-Sørensen similarity) among biomes using the total number of species. (C) Pairwise dissimilarity in trait hypervolumes (1-Sørensen similarity) among biomes using only those species that are considered as functionally similar and widespread. The lighter the cell the greater the dissimilarity.

analysis was beyond the scope of this study. More importantly, our results reinforce the importance of understanding how evolutionary and biogeographic history shape the functional diversity of biomes (Woodward et al., 2004; Pennington et al., 2006; Donoghue and Edwards, 2014; Moncrieff et al., 2016). The extensive sharing of species (and higher lineages) across biomes within, but largely not between, different biogeographic realms could serve both to homogenize functional diversity within realms and to provide clues about the characteristic traits that are selected for, or against, by different environments (Douma et al., 2012; Zanne et al., 2014, 2018).

Despite substantial hypervolume overlap among all the biomes (**Supplementary Figure 2**), tropical, temperate, and cold biomes all appear to occupy distinguishable regions of functional space (**Figures 5B,C**). The main traits differentiating biomes appear to be traits related to overall plant size, including both mature height and seed mass (**Supplementary Figures 7**–**10**), rather than leaf economics traits, as observed in more local, plot-based analyses (e.g., Douma et al., 2012). The exception is leaf P, which displayed substantial differences between tropical and temperate/polar biomes (**Supplementary Figures 7–10**), a pattern that has been observed in other analyses based on both species-specific values and whole ecosystem measurements (Kerkhoff et al., 2005; Swenson et al., 2012). Because leaf P was our most sparsely sampled trait, the differences between tropical and temperate/polar realms might be subject to biases from the imputation procedure. However, Leaf P does exhibit a significant phylogenetic signal (Kerkhoff et al., 2006), which suggests that the imputations should be unbiased (Swenson et al., 2017). Leaf P is also significantly higher in herbaceous growth forms (Kerkhoff et al., 2006), so the latitudinal shift from predominantly woody to herbaceous diversity (**Figure 4**) may also influence this pattern.

When we limited hypervolume analyses to only the most widespread and functionally common species in each biome, the individual biomes within each biogeographic realm overlapped less in functional trait space (compare **Figures 5B,C**), suggesting that these species may reflect phenotypes better adapted to particular environments. In this context, the xeric woodland biome is particularly interesting. In part due to the inclusion of the Chaco ecoregion, xeric woodlands cluster with the tropical biomes taxonomically (**Figure 5A**). But when the functional hypervolumes of all of the species are considered, they show much stronger similarity to the temperate biomes (**Figure 5B**). Finally, when only the most widespread and functionally common species are analyzed, xeric woodlands again show increased similarity to tropical biomes, but also they maintain high similarity with temperate grassland and Mediterranean biomes (**Figure 5C**). Transitions from warm, mesic environments to colder, drier, and more seasonal environments are facilitated by similar traits, e.g., herbaceous habit (Douma et al., 2012; Zanne et al., 2014, 2018), and like colder environments, communities in dry environments also tend to be more phylogenetically clustered (Qian and Sandel, 2017). Furthermore, the position of xeric zones at the boundary of the inter-tropical convergence zone makes them a geographical transition between the tropical and temperate realms. The fact that xeric woodlands are intermediate between the tropical and temperate realms both functionally and biogeographically further reinforces the idea that to better understand the functional diversity of biomes we must take into account their biogeographic and phylogenetic histories (Pennington et al., 2006; Donoghue and Edwards, 2014; Moncrieff et al., 2016).

## CONCLUSIONS

Any classification of terrestrial biomes imposes a small number of discrete categories on continuous gradients in climate and species distributions, and thus represents a gross simplification of the complex ecological landscape (Moncrieff et al., 2016, but see Silva de Miranda et al., 2018). Yet despite their potential for oversimplification, biomes are useful constructs for organizing and understanding the biodiversity and functioning of Earth's major terrestrial ecosystems, and trait-based approaches have high potential to help dynamically model global vegetation distributions.

In this study we have shown that much of the taxonomic diversity of all biomes represents species that are both narrowly distributed and functionally similar. Further, within biomes the most functionally distinctive species in each biome also tend to be geographically rare, while widespread species uniformly display low functional distinctiveness. Despite extensive taxonomic and functional overlap among biomes, they do cluster into distinguishable, biogeographically and climatically distinct units, especially when the functional clustering is based on the most widespread and functionally common species in each biome. However, advancing a functional understanding of biomes will require not just a better characterization of trait variation within and between vegetation types, but also information on the biogeographic and phylogenetic history of species assemblages and the relative abundance of species within biomes.

## DATA AVAILABILITY

Trait datasets and species range maps analyzed for this study can be download through the R package BIEN (Maitner et al., 2018). See http://bien.nceas.ucsb.edu/bien/ for more details about the BIEN database.

## AUTHOR CONTRIBUTIONS

SE-L, AK, BJE, DMN, and CV designed the study; SE-L and AK analyzed the data; SE-L and AK wrote the first version of the manuscript; CV, BB, NK, BSM, BM, RKP, BS, SS, J-CS, and SKW contributed to the collation and creation of the traits and range maps database. All authors contributed to the development and writing of the manuscript.

## FUNDING

SE-L, DMN, and AK were supported by a collaborative research grant from the US National Science Foundation (DEB-1556651). BJE was supported by National Science Foundation award DEB-1457812 and Macrosystems-1065861**.** CV was supported by the European Research Council (Grant ERC-StG-2014- 639706-CONSTRAINTS) and by the French Foundation for Research on Biodiversity (FRB; www.fondationbiodiversite.fr) in the context of the CESAB project Causes and consequences of functional rarity from local to global scales. J-CS considers this work a contribution to his VILLUM Investigator project Biodiversity Dynamics in a Changing World funded by VILLUM FONDEN (grant 16549). SKW was supported by the Strategic Science Investment Fund of the New Zealand Ministry of Business, Innovation and Employment's Science and Innovation Group.

## REFERENCES


## ACKNOWLEDGMENTS

We thank Dan Griffith and two reviewers for their constructive criticism which greatly improved this manuscript. We also thank all BIEN data contributors (see http://bien.nceas.ucsb.edu/bien/ data-contributors/ for a full list).

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2018.00219/full#supplementary-material

a means to quantitatively predict a broad range of species assemblages in NW Europe. Ecography 35, 364–373. doi: 10.1111/j.1600-0587.2011.07068.x


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Echeverría-Londoño, Enquist, Neves, Violle, Boyle, Kraft, Maitner, McGill, Peet, Sandel, Smith, Svenning, Wiser and Kerkhoff. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Patterns of Beta Diversity of Vascular Plants and Their Correspondence With Biome Boundaries Across North America

#### Jesús N. Pinto-Ledezma<sup>1</sup> \*, Daniel J. Larkin<sup>2</sup> and Jeannine Cavender-Bares <sup>1</sup> \*

<sup>1</sup> Department of Ecology, Evolution and Behavior, University of Minnesota, Saint Paul, MN, United States, <sup>2</sup> Department of Fisheries, Wildlife, and Conservation Biology, University of Minnesota, Saint Paul, MN, United States

Understanding why species composition and diversity varies spatially and with environmental variation is a long-standing theme in macroecological research. Numerous hypotheses have been generated to explain species and phylogenetic diversity gradients. Much less attention has been invested in explaining patterns of beta diversity. Biomes boundaries are thought to represent major shifts in abiotic variables accompanied by vegetation patterns and composition as a consequence of long-term interactions between the environment and the diversification and sorting of species. Using North American plant distribution data, phylogenetic information and three functional traits (SLA, seed mass, and plant height), we explicitly tested whether beta diversity is associated with biome boundaries and the extent to which two components of beta diversity—turnover and nestedness—for three dimensions of biodiversity (taxonomic, phylogenetic, and functional)—are associated with contrasting environments and linked to different patterns of historical climatic stability. We found that dimensions of vascular plant beta diversity are strongly coupled and vary considerably across North America, with turnover more influential in biomes with higher species richness and greater environmental stability and nestedness more influential in species-poor biomes characterized by high environmental variability. These results can be interpreted to indicate that in harsher climates with less stability explain beta diversity, while in warmer, wetter more stable climates, patterns of endemism associated with speciation processes, as well as local environmental sorting processes, contribute to beta diversity. Similar to prior studies, we conclude that patterns of similarity among communities and biomes reflects biogeographic legacies of how vascular plant diversity arose and was shaped by historical and ecological processes.

Keywords: beta diversity, biomes, functional traits, nestedness, niche conservatism, turnover, vascular plants

## INTRODUCTION

Ecological and evolutionary processes act in concert to modulate species traits that ultimately determine colonization and persistence of species within assemblages and across regions through evolutionary time (Simpson, 1953; Cavender-Bares et al., 2016). When species colonize novel areas and assemble into communities after speciation (Goldberg et al., 2011), their persistence depends

#### Edited by:

Colin Osborne, University of Sheffield, United Kingdom

#### Reviewed by:

Marta A. Jarzyna, The Ohio State University, United States Ana M. C. Santos, University of Alcalá, Spain

#### \*Correspondence:

Jesús N. Pinto-Ledezma jpintole@umn.edu Jeannine Cavender-Bares cavender@umn.edu

#### Specialty section:

This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution

Received: 14 May 2018 Accepted: 06 November 2018 Published: 27 November 2018

#### Citation:

Pinto-Ledezma JN, Larkin DJ and Cavender-Bares J (2018) Patterns of Beta Diversity of Vascular Plants and Their Correspondence With Biome Boundaries Across North America. Front. Ecol. Evol. 6:194. doi: 10.3389/fevo.2018.00194 on their capacity to deal with novel environmental conditions or track similar environments to which they are already preadapted (Reich et al., 2003; Emery et al., 2012; Cavender-Bares et al., 2016). Ultimately, the balance between ecological and evolutionary processes, in conjunction with physical (i.e., geological and climatic) events over long periods of time, determines patterns of taxonomic, functional, and phylogenetic biodiversity (Ricklefs, 1987; Emery et al., 2012; Cavender-Bares et al., 2016; Pinto-Ledezma et al., 2017; Símová et al., 2018), as well as the tendency for species to maintain their traits and environmental tolerances over evolutionary time ("niche conservatism"; Wiens and Donoghue, 2004; Crisp et al., 2009; Wiens et al., 2010).

A crucial challenge for understanding how biodiversity is structured and maintained through time is the integration of ecological and evolutionary processes at different spatial scales (Graham and Fine, 2008). To this end, analysis of beta diversity is a promising approach that allows heterogeneity in the distribution of gamma and alpha diversities to be quantified (Loreau, 2000; Baselga, 2010), thus permitting the evaluation of how species or lineages change across space and in response to environmental variation (Graham and Fine, 2008; Baselga, 2010). Here we ask how taxonomic, functional, and phylogenetic composition of plant species shift across biome boundaries and environmental gradients and consider factors that influence species' current distributions, including their evolutionary and biogeographic history. Although a number of studies have described large-scale patterns of beta diversity (variation in the composition of species) across ecological zones of vascular plants (e.g., Qian and Ricklefs, 2007; Qian, 2009) and vertebrates (e.g., Penone et al., 2016; Peixoto et al., 2017), very few studies explore species turnover (the replacement of some species by others) across zones (but see, McDonald et al., 2005). This study advances prior work by (1) deciphering large scale patterns of turnover and nestedness—as components of beta diversity for three dimensions of vascular plants diversity (taxonomic, phylogenetic, and functional) and (2) evaluating the roles of ecological and evolutionary processes in driving these observed patterns. In doing so, we use available databases of vascular plant species distributions and three critical functional traits (specific leaf area, seed mass, and plant height) related to resource acquisition, dispersal, and major axes of plant function to determine taxonomic, phylogenetic and functional distribution patterns of plants in North America. We then apply a partition procedure for beta diversity (Baselga, 2010; Leprieur et al., 2012) that allows the identification of species poor regions (differences in diversity between biomes—the nestedness component) and replacement regions (diversity interchange between biomes—the turnover component) (Baselga, 2010; Leprieur et al., 2012; see also **Box 1**).

Theories of niche conservatism posit that the colonization and adaptation of species to novel environments is uncommon. As consequence, most species accumulate within the environmental tolerances of their ancestors (Wiens and Donoghue, 2004; Wiens et al., 2010; Donoghue and Edwards, 2014), influencing assembly patterns of local communities (Ackerly, 2004; Cavender-Bares et al., 2009; Cavender-bares et al., 2018). In addition, even if corridors exist that connect areas with similar environmental conditions, species" persistence in newly colonized areas will depend on their functional traits and environmental tolerances, the result of legacies derived from evolutionary history within the context of past environments (Donoghue, 2008; Donoghue and Edwards, 2014; Cavender-Bares et al., 2016). Multiple lines of evidence demonstrate that niche conservatism is an important phenomenon across phylogenetic (Prinzing, 2001; Qian and Ricklefs, 2004) and spatial scales that limits the distributions of species (Crisp et al., 2009; Emery et al., 2012; Pinto-Ledezma et al., 2017).

Perhaps the most striking evidence for conservatism is conservatism of habitat preferences (Ackerly, 2004; Emery et al., 2012; Pinto-Ledezma et al., 2017) and biome associations within lineages (Crisp et al., 2009; Qian et al., 2017). Accordingly, evolutionary transitions among biomes are rare relative to speciation, and most lineages are restricted to the same biomes as their ancestors (Ackerly, 2004; Wiens and Donoghue, 2004; Crisp et al., 2009). However, an alternative view posits that biome shifts are common over evolutionary time, occurring mostly at biome transitions (Edwards and Donoghue, 2013; Donoghue and Edwards, 2014). When adaptation to new biomes has occurred, it may often have been the result of environmental changes that allow invasion of non-native species and subsequent in situ speciation (Donoghue, 2008; Donoghue and Edwards, 2014). Yet colonization and posterior adaptation to new biomes ultimately depend on legacy effects from trait and niche conservatism (Ackerly, 2004; Cavender-Bares et al., 2016). Finally, although lineage distributions might be affected by environmental changes, lineages expand or contract their ranges as a function of the concomitant expansion or contraction of the biomes with which they are associated (Crisp et al., 2009; Donoghue and Edwards, 2014; Cavender-Bares et al., 2016).

North America (United States and Canada) comprises 16.6% of the land surface of the Earth, spanning from 26◦ S to 85◦N and from 15◦W to 173◦E. This vast territory presents a wide diversity of climate regimes and land formations and encompasses many of the major ecosystem types and plant formations of the world (Graham, 1999; Qian, 1999). The major North American plant formations are considered to be the result of large-scale geological and climatic events during the Cenozoic (Braun, 1967; McKenna, 1983; Graham, 1999; Donoghue et al., 2001), as well as major biogeographic processes, including the colonization of temperate-adapted taxa from East Asia through the Beringia and tropical-temperate taxa from Europe via the North Atlantic land bridge (Graham, 1999; Donoghue et al., 2001). In terms of paleo-climatic changes, cooling and drying events in the Eocene (Graham, 2011) and glaciation cycles in the Pleistocene have important impacts on plant composition and diversity in North America (Davis, 1983; Clark et al., 1999).

These major plant formations, also known as "biomes" (Clements, 1936; Braun, 1967; Graham, 1999; Moncrieff et al., 2016), are regions with similar vegetation in terms of physiognomy, structure, and function and reflect the shared biogeographic history of past centers of diversification (Pennington et al., 2006; Donoghue and Edwards, 2014). Thus, from an evolutionary point of view, biomes represent large-scale

#### Box 1 | Beta diversity and its components

Beta diversity is a measure of the variation in species composition between two or more local communities or between local and regional communities (Koleff et al., 2003; Baselga, 2010) and can be partitioned and to understand the influence of species replacement or turnover (i.e., beta diversity caused by true substitution of species between sites) and richness differences or nestedness (i.e., beta diversity generated due to differences in species richness between areas) (Box 1: Figure 1) (Baselga, 2010; Legendre, 2014). Note that the same reasoning for variation in species composition between communities can be expanded to different dimensions of biodiversity (e.g., functional, phylogenetic) (Box 1: Figure 1). To illustrate both components, assume a given species pool (e.g., regional species richness) (Box 1: Figure 1). The species replacement or turnover component refers to the spatial change in species identities between two or more communities, and as turnover increases, focal communities tend to become more different from their neighbors in terms of composition, but maintain similar numbers of co-occurring species (Box 1: Figure 1; Equation 1B). Conversely, the richness difference or nestedness component refers to species losses or gains from one community to another, that is, nestedness happens when species-poor communities are strict subsets of species-rich communities (Box 1: Figure 1; Equation 1C). Although both turnover and nestedness can contribute separately to patterns of beta diversity, it is also possible for beta diversity to reflect the joint contributions of both components (Box 1: Figure 1).

This partitioning of beta diversity into turnover and nestedness allows the simultaneous evaluation of ecological and historical processes on current species diversity, as both components can result from different but not exclusive processes (Baselga, 2010; Leprieur et al., 2011, 2012; Legendre, 2014). For example, turnover results from historical events, geographical isolation, habitat specialization, and environmental filtering, while nestedness from selective forces (e.g., colonization, extinction) and environmental tolerances (Gaston et al., 2007; Baselga, 2010; Leprieur et al., 2011). One might therefore expect that species turnover occurs in regions that experienced less environmental variation over evolutionary time, as these regions have had more time for speciation, permitting the accumulation and maintenance of high numbers of species (Jansson, 2003; Wiens and Donoghue, 2004; Crisp et al., 2009). In addition, this environmental stability over time creates opportunities (e.g., through habitat connectivity) for species to colonize new areas with similar environmental conditions, increasing species turnover between regions (Stevens, 1989). In contrast, nestedness is expected to occur in regions that were severely affected by environmental changes in the past or regions that are subject to frequent disturbances (Baselga, 2010; Leprieur et al., 2011; Dobrovolski et al., 2012). Hence, beta diversity in these regions might result from past extinction events and the recent colonization of species that are capable of tolerating harsh environmental conditions (Jansson, 2003; Baselga, 2010; Legendre, 2014).

The general equations that we used are as follows:

Equation, 1A Beta diversity based on Sørensen dissimilarity index (βsor)

$$
\beta\_{\text{SOI}} = \frac{b+c}{(2\mathbf{a} + \mathbf{b} + c)}
$$

Equation, 1B Turnover component (βsim)

$$
\beta\_{\text{crit}} = \frac{\min(b, c)}{\mathfrak{a} + \min(b, c)}
$$

Equation, 1C Nestedness component (βnes)

$$\beta\_{\text{tree}} = \frac{\mathsf{m}\mathsf{new}(\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\beta}}\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\beta}}\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\mathsf{\prime$$

where a is the number of species in both communities, b is the number of species exclusive to a focal community and c is the number of species exclusive to the adjacent communities. These equations were used for taxonomic beta diversity calculations. The equations used for the calculation of phylogenetic and functional beta diversity can be found in the Supplementary Material.

Figure 1 | Conceptual framework for a partition procedure of beta diversity into nestedness and turnover components. Turnover is the difference between two assemblages due to contrasting subsets of species or lineages from a source pool, potentially caused by environmental sorting processes or speciation processes associated with endemism forced by topographic barriers, while nestedness is the difference in species composition between two assemblages due to attrition of species in one assemblage relative to the other. The example involves the three dimensions considered in our study and three biomes. Also, as species carry both functional and phylogenetic information (Safi et al., 2011; Swenson, 2011), it is likely that observed patterns of nestedness and turnover will be equivalent for the three dimensions. (A) All biomes have the same number of species and display a pattern of species turnover within and between biomes. (B) Shows complete nestedness (differences in diversity between biomes), where biomes 2 to 3 represent a subset of more diverse biomes (biome 1), potentially as a response to environmental filters. (C) Shows the joint contribution of both components to the patterns of beta diversity.

ecological niches for most plant lineages (Pennington et al., 2006; Crisp et al., 2009; Donoghue and Edwards, 2014). Although biomes are not intended to correspond precisely with climatic regions, their distributions are strongly tied to climate and other environmental factors (Moncrieff et al., 2016) and their boundaries arise due to differences in abiotic conditions (i.e., geological, climatic, and edaphic) and biotic similarity (i.e., compositional and functional dimensions of diversity) between biomes (Donoghue and Edwards, 2014; Moncrieff et al., 2016).

Given that ecological and historical processes can either increase or decrease co-occurrence and phylogenetic relatedness within communities and between biomes (Ricklefs, 1987; Cavender-Bares et al., 2009; Crisp et al., 2009), alternative explanations can be used to describe patterns of beta diversity for vascular plants in North America or other regions. According to the evidence for biome conservatism (Crisp et al., 2009; Qian et al., 2017), we expect differential effects of ecological and historical processes among biomes, resulting in differential contributions of beta diversity components, i.e., turnover and nestedness, on assembly of communities within and across biomes (see **Box 1**). Consider that most vascular plants originated during the Ordovician period (∼470 Ma), when the environment was hot and wet; cold and drought tolerances likely evolved later during the Eocene (∼50 Ma) after the start of the global cooling period (Zachos et al., 2008; Zanne et al., 2014; Cavender-Bares et al., 2016). We therefore predict that turnover in species should be the dominant contributor to beta diversity in biomes where the environment is most stable, as environmental stability over time allows the creation of large species pools and consequently the interchange of species between biomes (**Box 1**) and, conversely, that nestedness (**Box 1**) should be most important in less stable environments, such as cold and dry environments most impacted by glaciation cycles and associated climatic changes, which could have led to loss of species (Qian et al., 2005; Qian and Ricklefs, 2007; Graham and Fine, 2008; Dobrovolski et al., 2012; **Figure S1**). We further expect functional beta diversity to follow the same predicted pattern as taxonomic and phylogenetic diversity (Safi et al., 2011; Swenson, 2011; Penone et al., 2016) if variation in functional traits that underlie the assembly of communities is associated with environmental variation and if those traits are phylogenetically conserved (Cavender-Bares et al., 2009, 2016).

In addition, given that functional traits facilitate dispersal, establishment, and persistence of plants across environmental gradients (Westoby, 1998; Garnier et al., 2001; Moles et al., 2005; Qian, 2009; Zanne et al., 2014; Díaz et al., 2015; Moles and Gibson, 2017; see also **Figure S2**), these traits can modulate patterns of beta diversity within communities and between biomes (Qian, 2009; Cavender-Bares et al., 2016). Accordingly, it is expected that patterns in functional turnover and nestedness should show strong relationships with environmental gradients (Qian et al., 2005; Gaston et al., 2007; Qian and Ricklefs, 2007; Qian, 2009) and consequently with biome distribution (Díaz et al., 2015; Moles and Gibson, 2017). Therefore, we expect the influence of trait turnover and nestedness to vary between biomes with respect to traits related to dispersal, establishment, and persistence. For example, cold and dry biomes dominated by species with high dispersal ability (e.g., small seeds), high persistence capacity (e.g., high specific leaf area), and high resistance to environmental disturbances (e.g., low plant height or non-woody growth forms) (**Figure S2**) might be dominated by the nestedness component (Qian, 2009; Dobrovolski et al., 2012), given that vagile species are less affected by environmental barriers and non-woody species are better able to colonize and persist in areas subject to harsh environmental conditions (Garnier et al., 2001; Qian, 2009; Zanne et al., 2014; Díaz et al., 2015). Conversely, less seasonal biomes (hotter and wetter) may be dominated by less vagile species with low persistence capacity and low resistance to environmental change (**Figure S2**) and show greater turnover given that these biomes were less influenced by drastic environmental changes over time (**Box 1**). These biomes would thus have allowed greater accumulation of species within them, consequently producing larger species pools for assembly into local communities (Crisp et al., 2009; Cavender-Bares et al., 2016).

## MATERIALS AND METHODS

## Distributional, Phylogenetic, and Functional Data

We obtained species distribution data for vascular plants from two sources. Species ranges were downloaded from the BIEN database (Botanical Information and Ecology Network, http://bien.nceas.ucsb.edu/bien/, last accessed February 21, 2018) (Enquist et al., 2016) and occurrence points were downloaded from the GBIF database (Global Biodiversity Information Facility, https://www.gbif.org, last accessed February 27, 2018). Using the occurrence points we built species ranges for species with at least 10, non-disjunct occurrence points (Davis-Rabosky et al., 2016; Meyer et al., 2018). We did not consider species with <10 points because ranges constructed with too few occurrences tend to overestimate range size (Meyer et al., 2018). The phylogenetic hypothesis used was obtained from the recently published Spermatophyta mega-phylogeny (Smith and Brown, 2018). This phylogeny was reconstructed using two different backbones [i.e., OTB (Open Tree of Life backbone) and MB (Magallón backbone)], and missing species were added according the Open Tree of Life (see Smith and Brown, 2018 for details). Although in this phylogeny many taxa are unresolved, it is the most comprehensive phylogeny for plant species to date, containing a total of 353,185 and 356,305 species for the OTB and MB backbones, respectively (Smith and Brown, 2018). In this study, we used the mega-phylogeny constructed under the MB backbone. Functional trait data were obtained from the TRY (Kattge et al., 2011) and BIEN (Enquist et al., 2016) databases and comprised seed mass (SM), whole plant height (WPH), and specific leaf area (SLA). These three traits were selected because they are related to different processes, including dispersal, establishment, and persistence (Westoby, 1998; Díaz et al., 2015; Moles and Gibson, 2017). We included only species that had all three kinds of data (occurrences, phylogenetic information, and traits). This yielded a dataset with 12,317 species, ∼80% of the nearly 15,500 vascular plants registered for North America (Ulloa-Ulloa et al., 2017).

From these data, we rasterized vascular plants' geographical ranges into a grid of 1◦ × 1 ◦ resolution (approximately 110 × 110 km near the equator), maintaining cells with more than 25% of land area to obtain a species presence-absence matrix (PAM) consisting of assemblages within grid cells (as rows) × species (as columns). We used this resolution because it is commonly used in macroecological studies of beta diversity (e.g., McKnight et al., 2007; Melo et al., 2009; Dobrovolski et al., 2012) and is appropriately scaled to account for errors related to species ranges (Kreft and Jetz, 2007). The number of species per cell (assemblage) were obtained by summing all columns; assemblages with fewer than 10 species were removed due to their potential to produce spurious results (Kreft and Jetz, 2007). The resulting PAM included a total 3,298 assemblages and 12,317 species. Data were retrieved using the R packages BIEN (Maitner et al., 2018) and spocc (Scott et al., 2017), and the PAM was built using the R package letsR (Vilela and Villalobos, 2015).

## Environmental Data

We obtained environmental data for annual mean temperature (BIO1), maximum temperature of the warmest month (BIO5), minimum temperature of the coldest month (BIO6), annual precipitation (BIO12), precipitation in the driest quarter (BIO9), precipitation in the warmest quarter (BIO18), and altitude (Alt) from the 10 arc-min WorldClim database (Hijmans et al., 2005). We also obtained two derived climate indices: actual evapotranspiration (AET) and potential evapotranspiration (PET). These indices reflect fluxes of water related to habitat productivity and the proportional drying power of the environment (Qian et al., 2005). These environmental measures are standard variables for describing climatic characteristics of species ranges and are commonly used in macroecological studies of plants (Qian et al., 2005; Qian and Ricklefs, 2007; Moles et al., 2014). To reduce collinearity and variability of environmental variables we used principal component analysis (PCA). Prior to implementing PCA, we standardized all variables to a mean of 0 and standard deviation of 1. We then performed PCA and extracted mean values for each assemblage for each PC axis using the R package raster (Hijmans et al., 2017).

## Beta Diversity Measures

Beta diversity was measured for each dimension of vascular plant diversity (i.e., taxonomic, phylogenetic, and functional) using the Sørensen dissimilarity index (Sørensen, 1948), and a partitioning procedure that takes numbers of species into account (Baselga, 2010; Leprieur et al., 2012). The equations used to calculate taxonomic and phylogenetic beta diversity and their components are detailed in Baselga (2010) and Leprieur et al. (2012), respectively and in the (**Supplementary Data 1**). Briefly, phylogenetic diversity is calculated as the total branch length of a tree that involves all co-occurring species in a community, and functional beta diversity uses the same formulae using functional dendrograms for the same set of species (see below). The partitioning procedure (**Box 1**) separates beta diversity (βsor, **Box 1**: Equation 1A) into two components: the turnover component (βsim, **Box 1**: Equation 1B) and the nestedness component (βnes, **Box 1**: Equation 1C), where βsim represents the true spatial species replacement between assemblages, without the influence of richness gradients and βnes represents the difference in species richness between assemblages that occurs when species-poor assemblages are subsets of species-rich assemblages (McKnight et al., 2007; Baselga, 2010; Leprieur et al., 2012).

Prior to calculating functional beta diversity using the three functional traits (i.e., SM, WPH, SLA) we constructed a functional dendrogram that estimates species functional differences or dispersion in functional space (Petchey and Gaston, 2006; Safi et al., 2011). The functional dendrogram used Euclidean pairwise distances and the Unweighted Pair Group method with Arithmetic Mean (UPGMA) clustering algorithm (Petchey and Gaston, 2006; Safi et al., 2011). We chose Euclidean distances because our functional traits are continuous data. The resulting functional dendrogram represents a continuous measure of functional diversity (FD), where high FD indicates high species differences in functional space, whereas low FD indicates that species are more similar in functional space (Petchey and Gaston, 2006). FD was determined by summing the branch lengths of the functional dendrogram that connects all species co-occurring in the community (Safi et al., 2011). Subsequently, functional beta diversity was calculated using the functional dendrogram. The same equations used to calculate the phylogenetic beta diversity (see Equations 2A–C in the **Supplementary Material**) were used for functional beta diversity.

Similar to the multivariate trait dendrogram, we built individual trait dendrograms and repeated all functional beta diversity calculations for each individual trait (i.e., SM, SLA, and WPH), given that these traits do not necessarily covary and that each may be particularly relevant to a specific function. For example, seed mass (SM) influences dispersal, with small-seeded plants generally having higher dispersal capacity than large-seeded plants; however, larger seeds are typically better provisioned and may confer advantages in establishment (Westoby, 1998; Moles et al., 2005; Qian, 2009). Specific leaf area (SLA) is related to resource acquisition (Reich, 2014), photosynthetic capacity (Wright et al., 2004), and growth and competitive ability and is positively correlated with relative growth rate (RGR) across species (Garnier et al., 2001). Finally, whole plant height (WPH) is associated with establishment and resistance to environmental disturbances (Moles et al., 2009; Moles and Gibson, 2017).

We calculated beta diversity between pairs of assemblages (1◦ × 1 ◦ cells) using a moving window approach (Melo et al., 2009; Dobrovolski et al., 2012). A focal assemblage and its neighbors were selected using defined window sizes and beta diversity was calculated as the mean beta diversity between each focal assemblage and its neighbors (see **Figure S3**). To account for spatial dependency in the analysis we calculated beta diversity using window sizes varying from 1◦ to 5◦ degrees in 1◦ steps. For each window size, we calculated taxonomic, phylogenetic, and functional beta diversity (hereafter βsor−tax, βsor−phy, βsor−func). The different window sizes produced similar patterns of beta diversity (**Figure S3**), and we present our results using a 2◦ window size. All calculations were performed in R v3.4 (R Core Team, 2018) using customized scripts and core functions from the packages betapart (Baselga and Orme, 2012), and CommEcol (Melo, 2017).

Using the beta diversity calculations (i.e., βsor−tax, βsor−phy, βsor−func), we calculated the ratio between βnes and βsor (hereafter "βratio") and the deviations of BDphy given the BDtax (hereafter "βdev"). The βratio for the three dimensions of vascular plants were calculated separately as βratio = βnes βsor and identifies the relative influence of the βsim or βnes components of beta diversity (Dobrovolski et al., 2012; Peixoto et al., 2017). Values < 0.5 indicate that beta diversity is determined mostly by species replacement between assemblages (i.e., βsim or turnover component); conversely, values > 0.5 indicate that beta diversity is mostly influenced by the nestedness component (βnes). βdev were calculated as βdev = βsor−tax − βsor−phy βsor−tax and reflects the degree of lineage exchanges between areas over time (Graham and Fine, 2008; Peixoto et al., 2017), where high positive values suggest little lineage exchange (βsor−tax > βsor−phy) and negative values suggest high lineage exchange (βsor−phy > βsor−tax) (Peixoto et al., 2017).

## Statistical Analyses

We evaluated the structure and spatial distribution of beta diversity using spatial correlograms (Diniz-Filho et al., 2003) and spatial congruence of the different dimensions of diversity (taxonomic, phylogenetic, and functional) using linear correlations. We controlled for spatial autocorrelation using Clifford's method to obtain effective degrees of freedom for Pearson's coefficients (Clifford et al., 1989; Pinto-Ledezma et al., 2017). We evaluated the influence of the environmental variables on beta diversity using ordinary least-square (OLS) models. Because regression residuals tend to display high autocorrelation in spatially structured data (Diniz-Filho et al., 2003), we also modeled the relationship between beta diversity and the environmental predictors using simultaneous autoregressive models (SAR) (Kissling and Carl, 2008). SAR models account for spatial autocorrelation by adding an extra term (autoregressive) in the form of a spatial-weight matrix that specifies the neighborhood of each assemblage and the relative weight of each neighbor (Kissling and Carl, 2008). Both OLS and SAR models estimate single coefficients accounting for cells in a geographic domain. In addition, to explore the local influence of the environmental variables on beta diversity we used geographically weighted regressions (GWR) (Fotheringham et al., 2002). All statistical analyses were performed in R using the packages sp (Pebesma and Bivand, 2005), spdep (Bivand and Piras, 2015), spgwr (Bivand et al., 2017), and SpatialPack (Vallejos and Osorio, 2014).

## RESULTS

Beta diversity (βsor) and its components (βnes and βsim) were unevenly distributed across North America (**Figure S3**). βsim for all three dimensions of vascular plant diversity was higher in dry and temperate biomes in southern North America and in polar biomes in northern North America, whereas, βnes was higher at mid- to high latitudes positioned between polar and temperate biomes (**Figure 2**). Spatial structure of the βnes and βsim components were similar across the three dimensions of plant vascular diversity (**Figures 2A–C**; **Figure S3**), which showed high spatial congruence based on spatial correlations (**Table 1**) and correlograms (**Figure S4**). Interestingly, patterns of turnover and nestedness for specific functional traits (**Figures 2D–F**) were more divergent than were taxonomic, phylogenetic, and multivariate functional axes of diversity. For example, βsim for whole plant height was higher within and between deciduous forest and temperate savanna (**Figure 2D**), while βsim for seed mass was higher within the transitions of deciduous forest, savanna, and scrubland (**Figure 2E**). For specific leaf area, the pattern was even more complex, with higher βsim within and between temperate grasslands and deserts at the west and within and between temperate forest and grasslands at the east (**Figure 2F**). In addition, we also identify high βsim between the coniferous forest biome and tundra biome for the three dimensions (**Figures 2A–C**) and for seed mass and specific leaf area (**Figures 2E,F**).

These results support the hypothesis that historical processes (e.g., speciation) were critical in driving beta diversity in North American vascular plants (**Figure 2**). Following our predictions, we found that βnes (differences in species richness between biomes) was more important at mid- to high-latitudes in North America (blue areas in **Figure 2A–C**), a pattern consistent with the expectation that biomes at these latitudes (∼50◦ -60◦ N) were more recently colonized. Conversely, <sup>β</sup>sim (true species substitution between biomes) contributed more strongly to species and lineage composition in southern North American biomes (red areas in **Figures 2A–C**). The pattern is consistent with the interpretation that much of the vascular plant species diversity at these latitudes (below 40◦ N) has persisted and accumulated over evolutionary time scales, resulting in higher concentrations of species (**Figure S1C**). For example, Mediterranean and deciduous forest biomes attain up to ∼4,000 and ∼3,200 species, respectively. The finding of high βnes in vascular plants at northern latitudes (arctic biomes), which is consistent with the interpretation of more recent colonization, is corroborated by the pattern of low βdev (dark blue cells in **Figure S5**), indicating low lineage exchange, such that relatively few lineages were able to colonize and subsist within these regions.

Finally, we found that environmental conditions had differential associations with patterns of nestedness and turnover (**Table 2**). βnes for taxonomic, phylogenetic, and multivariate functional dimensions of plant biodiversity was negatively associated with PC1 (related to variation in temperature; see **Table S1** for a summary of the PCA) and positively associated with PC2 (related to variation in precipitation), PC3 (related to variation in seasonality), and altitude, although variation in βnes was mainly explained by PC3 (**Table 2**). βsim was more strongly associated with PC2 and PC1. These results are most evident when the local influence of the environment on beta diversity components are examined (**Figure 3**). GWR models indicated that the environment had strong effects on beta diversity components for all three dimensions of diversity

[Quasi-global R 2 for βnes (Taxonomic = 0.836; Functional = 0.875; Phylogenetic = 0.880) and βsim [Taxonomic = 0.971; Functional = 0.883; Phylogenetic = 0.908)], although the environmental influence over βnes and βsim varied within and

βnes (blue color cells) and values < 0.5 strong influence of βsim (red color cells).

## DISCUSSION

Patterns of species diversity and composition within local communities and among regions are shaped by different historical and ecological processes (Ricklefs, 1987; Chesson, 2000; Cavender-Bares et al., 2009). In this study, we used large datasets of vascular plant diversity in North America and applied procedures partitioning beta diversity (Baselga, 2010; Leprieur et al., 2012) to demonstrate the simultaneous influence of historical and ecological processes on large-scale patterns of species turnover and nestedness.

More specifically, we found high spatial correlation among taxonomic, functional, and phylogenetic dimensions of diversity, suggesting that these axes are coupled and may be driven by similar processes (**Figures 2A–C**). We also found that vascular plant beta diversity varies considerably across North America, with turnover more influential in biomes with higher species richness and greater environmental stability and nestedness

between biomes across North America (**Figure 3**).

more influential in species-poor biomes characterized by high environmental variability (**Figures 2**, **3**; **Figures S1**, **S3**). In other words, large-scale, long-term evolutionary and environmental processes are reflected in contemporary assembly of local and regional floras across North America. Further, when functional traits are analyzed separately, new insights regarding the assembly of local communities and the exchange of species and lineages among biomes emerge. These findings are consistent with other studies addressing taxonomic diversity of vascular plants (Qian et al., 2005; Qian and Ricklefs, 2007; Qian, 2009), phylogenetic diversity of angiosperms (Qian et al., 2013), functional diversity of trees (Siefert et al., 2013), and multiple dimensions of vertebrate diversity (McDonald et al., 2005; McKnight et al., 2007; Melo et al., 2009; Dobrovolski et al., 2012; Penone et al., 2016; Peixoto et al., 2017).

As expected, our results showed differential contributions of turnover and nestedness to beta diversity of vascular plant communities both within and between biomes (**Figure 2**). The differential contributions of these components can be explained by alternative but non-mutually exclusive processes.

TABLE 1 | Spatial correlations for beta diversity and its components (turnover and nestedness) for the three dimensions of biodiversity (taxonomic, functional, and phylogenetic).


Degrees of freedom (d.f.) were corrected using Clifford's method (Clifford et al., 1989). r, Pearson's correlation coefficient; P, associated p-value.

For instance, biome niche conservatism (BNC) predicts that evolutionary biome transitions are infrequent and that most lineages remain in the niche conditions of their ancestors (Crisp et al., 2009; Qian et al., 2017). Accordingly, our results showed a shift from turnover (true species substitution between biomes) in warmer and wetter biomes (**Figure 3**; **Figure S1D**) to nestedness (differences in species richness between biomes) in colder and drier biomes (**Figures 2**, **figures 3A–C**, **Figure S1D**). While conservatism of species climatic tolerances is one likely explanation for these patterns, it also possible that geographical barriers may have caused the observed nestedness pattern; montane areas encompass high environmental variability within short distances and may thus constrain species distributions (Melo et al., 2009; Dobrovolski et al., 2012). These findings are consistent with those of Qian and Ricklefs (2007) who showed that turnover in vascular plants decreases as latitude increases (see also Qian et al., 2005).

This change in the contribution of turnover and nestedness could imply that most vascular plants remained—and thus had more time for speciation—within their ancestral biome conditions (warmer and wetter) (**Figures 2**; **Figure S1C**; Crisp et al., 2009; Hawkins et al., 2011; Qian et al., 2017). In fact, most vascular plant families originated and diversified during warmer/wetter conditions in the Ordovician (∼470 Ma) and Cretaceous (∼140 Ma) periods (Magallón et al., 2015; Vamosi et al., 2018), which potentially allowed the accumulation of more species within tropical-like biomes (Fine and Ree, 2006; Hawkins et al., 2011) and consequently contributed to the formation of larger species pools (Chesson, 2000; Crisp et al., 2009; Cavender-Bares et al., 2016). Further, during most of its evolutionary history, North America was colonized by plant lineages from eastern Asia and western Europe via the Beringia and North Atlantic land bridges, respectively (McKenna, 1983; Graham, 1999; Donoghue et al., 2001). These land bridges allowed some lineages (e.g., Castanopsis and Lithocarpus within the family Fagaceae) to track similar climatic tolerances in different continents, as their climatic tolerances expanded over evolutionary time (Donoghue, 2008; Cavender-Bares et al.,


\*P < 0.05; \*\*P < 0.01; \*\*\*P < 0.001; NS, non-significant. PC1 was correlated with temperature variables, PC2 with precipitation variables, and PC3 with additional variation in temperature and precipitation.

2016). The rapid turnover within these biomes (**Figures 3A–C**, **Figure S1**) could also have emerged due to the high habitat specialization enabled by less-seasonal biomes (Stevens, 1989; Qian, 1998; Qian and Ricklefs, 2007), which promotes high species co-occurrence through the reduction of interspecific competition (Usinowicz et al., 2017).

In contrast to the warmer/wetter biomes, colder/drier biomes (e.g., coniferous forest, tundra) had lower species co-occurrence (**Figure S1C**) and were characterized by higher proportions of gymnosperms (e.g., conifers), shrubs, and herbaceous species (Qian, 1998; Graham, 1999). Our results show strong shifts across all dimensions of diversity from turnover to nestedness (**Figures 2A–C**) concurrent with the transition from southern North American biomes (i.e., deciduous forest, grasslands, and shrublands and deserts) to the coniferous forest biome in the north (**Figures 2A–C**). This change from turnover to nestedness aligns with a shift from species-rich to species-poor biomes (**Figure S1C**; see also Qian, 1998), which is likely driven by recent glaciation in North America and postglacial colonization by species preadapted to severe environmental conditions (Garnier et al., 2001; Qian and Ricklefs, 2007; Crisp and Cook, 2011). Indeed, most species in the tundra biome comprise preadapted subsets of the species pool from adjacent biomes (Graham, 1999). Another possible explanation is that this pattern would emerge from evolutionary shifts in species' functional traits related to climatic tolerances that allowed the colonization of and in situ speciation of preadapted taxa (e.g., conifers) during cooling periods in the Oligocene-Miocene (Crisp and Cook, 2011). These findings suggest that long-term environmental change and harsh environmental conditions play important roles in filtering species pools within these biomes (Graham and Fine, 2008; Cavender-Bares et al., 2009; Cavender-bares et al., 2018).

We were also interested in exploring patterns of beta diversity of specific functional traits related to resource acquisition, dispersal, and major axes of plant function (**Figure S2**). We hypothesized that different traits would show differential sensitivity to environmental conditions within and between biomes. Accordingly, spatial pattern of turnover and nestedness components in individual traits (**Figures 2D–F**) did vary from the pattern exhibited for multivariate functional beta diversity (**Figure 2C**), though both exhibited decrease in turnover as environmental conditions become colder and drier (**Figures 2D–F**, **Figure S1C**). This suggests that environmental filtering led to convergent functional composition in communities that shared similar environmental conditions (Cavender-Bares et al., 2009; Swenson et al., 2012; Siefert et al., 2013; Símová et al., 2018). Moreover, several lines of evidence indicate that plant functional traits are correlated with environmental conditions (Swenson and Weiser, 2010; Swenson et al., 2012; Moles et al., 2014; Símová et al., 2018).

We also expected the relative contribution of turnover and nestedness to beta diversity in functional trait composition to differ across biomes. Our results showed that trait dissimilarity for the three functional traits considered in this study across biomes was driven mainly by differences in functional richness/composition (nestedness) rather than functional substitution (turnover) (**Figures 2D–F**). These patterns were most evident in biomes characterized by harsh environmental conditions and low species richness (**Figure S1**). Studies of vascular plants (Qian and Ricklefs, 2007; Qian, 2009) and vertebrates (Dobrovolski et al., 2012; Peixoto et al., 2017) indicate that recent, post-glacial range expansions toward northern latitudes have been constrained by species' climatic tolerances and by traits related to dispersal and persistence capacity (Qian and Ricklefs, 2007; Qian, 2009; Dobrovolski et al., 2012, but see Penone et al., 2016 for an alternative view). In general, our results are congruent with these studies; we found that northern tundra and coniferous biomes were dominated by species with high vagile capacity, high persistence capacity, and high resistance to environmental disturbances (**Figure S2**) (Moles et al., 2006, 2011; Díaz et al., 2015; Símová et al., 2018.

Our analyses also showed high functional turnover for plant height (related to resistance to environmental disturbances) and specific leaf area (related to persistence capacity) within and between biomes of temperate deciduous forest, temperate grasslands, and deserts (**Figures 3D,F**). Of particular interest is that these patterns of functional turnover appear to have been driven not by numbers of co-occurring species (**Figure S1C**), but by the evolutionary legacy of functional traits that enabled interchange of species between biomes with similar environmental conditions over time (Crisp et al., 2009; Cavender-Bares et al., 2016). This observation is based on the hypothesis that similar environmental conditions act as selective forces that filter species with similar functional traits (Penone et al., 2016) and the spatial distribution of functional traits across biomes (**Figure S2**). For example, biomes characterized by warmer climates are composed mostly by species with lower values of SLA (**Figure S2C**).

While considerable attention has been paid to understanding patterns of species diversity, much less emphasis has been placed on understanding patterns of beta diversity and its constituent components. Our findings with respect to vascular plant diversity in North America show that partitioning of beta diversity into turnover and nestedness components reveals how similarity among communities and biomes likely reflects historical processes and biogeographic legacies that have influenced how vascular plant diversity arose and assembled ecologically. Specifically, our results suggest that turnover is mainly controlled by tolerances environmental conditions over evolutionary time (e.g., limited niche evolution) that constraint the dispersal and adaptation of different niche conditions away from their ancestors, while nestedness seems an outcome of recent colonization events that occurred after last glaciation events during the Pleistocene. Finally, although it has been suggested that biome boundaries represent opportunities for species interchange (Moncrieff et al., 2016), biome boundaries are not static as their distribution is tightly linked to dynamic climatic and geological processes (Braun, 1967; Williams et al., 2004; Williams and Jackson, 2007; Cavender-Bares et al., 2016). Thus, it is important to consider biomes as dynamic entities that have changed in terms of distribution, structure and composition over time (Williams et al., 2004) and that future climate change will alter distribution of vascular plant diversity and thus the current biome distributions and boundaries (Williams et al., 2004; Williams and Jackson, 2007).

## AUTHOR CONTRIBUTIONS

JC-B, JP-L, and DL conceived the study and planned the analyses; JP-L managed the project. JP-L compiled all data used in the paper and conducted the analyses. JP-L led the writing and all authors contributed throughout the whole process.

## FUNDING

JP-L was supported by the University of Minnesota College of Biological Sciences' Grand Challenges in Biology Postdoctoral Program. NSF/NASA Dimensions of Biodiversity program (DEB-1342872).

## ACKNOWLEDGMENTS

We thank Colin Osborne, Daniel M. Griffith, and Christopher Still for the invitation to submit a study to the special issue Revisiting The Biome Concept With A Functional Lens and Fabricio Villalobos for his advice with beta diversity analyses.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2018.00194/full#supplementary-material

## REFERENCES


open habitats in driving species richness gradients. J. Biogeogr. 44, 1683–1693. doi: 10.1111/jbi.12939


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Pinto-Ledezma, Larkin and Cavender-Bares. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Molecular Phylogeny and Biogeography of Adenocaulon Highlight the Biogeographic Links between New World and Old World

Tao Deng1†, Yongsheng Chen2†, Hengchang Wang<sup>3</sup> , Xiaoshuang Zhang1, 4, Sergei Volis <sup>1</sup> , Ziyoviddin Yusupov <sup>1</sup> , Hong Qian<sup>5</sup> and Hang Sun<sup>1</sup> \*

*<sup>1</sup> Key Laboratory for Plant Diversity and Biogeography of East Asia, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming, China, <sup>2</sup> College of Life Science, Northeast Agricultural University, Harbin, China, <sup>3</sup> Key Laboratory of Plant Germplasm Enhancement and Specialty Agriculture, Wuhan Botanical Garden, Chinese Academy of Sciences, Wuhan, China, <sup>4</sup> Key Laboratory of Plant Resources Conservation and Utilization, College of Biology and Environmental Sciences, Jishou University, Jishou, China, <sup>5</sup> Research and Collections Center, Illinois State Museum, Springfield, IL, United States*

#### Edited by:

*Daniel M. Griffith, Oregon State University, United States*

#### Reviewed by:

*Tao Su, Xishuangbanna Tropical Botanical Garden (CAS), China Xiaolei Huang, Fujian Agriculture and Forestry University, China*

\*Correspondence:

*Hang Sun sunhang@mail.kib.ac.cn*

*† These authors have contributed equally to this work.*

#### Specialty section:

*This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution*

Received: *26 August 2017* Accepted: *29 November 2017* Published: *04 January 2018*

#### Citation:

*Deng T, Chen Y, Wang H, Zhang X, Volis S, Yusupov Z, Qian H and Sun H (2018) Molecular Phylogeny and Biogeography of Adenocaulon Highlight the Biogeographic Links between New World and Old World. Front. Ecol. Evol. 5:162. doi: 10.3389/fevo.2017.00162* *Adenocaulon* (Asteraceae) is a small genus with only five species but has a broad amphi-Pacific distribution pattern with three species distributed disjunctly in South America, Central America, and North America and two endemic species spanning from eastern Asia to the Himalayas. To trace the biogeographic pattern of the genus, we reconstructed its phylogenetic relationships and diversification history based on one nuclear and eight plastid gene regions. Our results showed that *Adenocaulon* is monophyletic and may have originated in Central America during the Miocene, dispersed into North America and finally reached the Himalayas via the Bering Land Bridge. The hypothesized trajectory implies that long-distance dispersal may have played an important role in the formation of the distribution of this group of species. This hypothesis seems to have gained support from the special morphological structure of fruits of the genus.

Keywords: Adenocaulon, Asteraceae, East Asia, America, herbaceous flora

## INTRODUCTION

The genus Adenocaulon Hook., consisting of five species, is a small genus in the angiosperm family Asteraceae, and has a broad amphi-Pacific distribution (Bittmann, 1990). Two species are endemic to East Asia while the other three occur in North and South America. Among them, Adenocaulon nepalense Bittmann is native to Bhutan, Nepal and India, and Adenocaulon himalaicum Edgew is distributed from the Himalayas to northeastern Asia and Japan (Bittmann, 1990; Funk and Hind, 2016). In the New World, Adenocaulon bicolor Hook. is distributed in North America, Adenocaulon lyratum S. F. Blake in Guatemala and Chiapas of Mexico, and Adenocaulon chilense Less. in Chile and Argentina (**Figure 1**; Funk and Hind, 2016). Thus the five species generally present a disjunct distribution pattern, with no species of Adenocaulon present in Africa, the Mediterranean, tropical Asia or Oceania.

Historically the taxonomic position of Adenocaulon has undergone extensive circumscription. Jensen and Kim (1996) placed it in the tribe Mutisieae, Katinas (2000) placed it in Cichorioideae, and Kim et al. (2002) transferred the genus to the Nassauvieae. A recent molecular phylogenetic study strongly supports the two American species, A. bicolor and A. chilense, forming a single lineage within Mutisieae (Panero and Funk, 2008), with the relationship of Adenocaulon to the

other genera in Mutisieae being unresolved. The phylogenetic relationship and geographic pattern of the genus are still unclear. The tribe Mutisieae, as well as other members of the early diverging branches of Asteraceae (e.g., Barnadesioideae and Mutisioideae s.l.), has a predominant center of diversity in South America, which is believed to be the most possible ancestral area of these lineages (Katinas et al., 2009; Ortiz et al., 2009).

A distinct feature of members of Asteraceae is one-seeded cypselae fruit. Further, the majority of extant Asteraceae species have a modified calyx above the inferior ovary structure called a pappus, allowing for seed dispersal by wind and biological agents over great distances (Katinas et al., 2013). However, this specialized morphological trait is not evident in Adenocaulon. Members of this genus usually grow in under-stories of temperate moist forests (Katinas, 2000), which may well be a good case study of the biogeographic link between the Old World (specifically East Asia) and the New World by inferring the dispersal mechanism and biogeographic patterns across the Pan-Pacific areas.

Although the East Asian and North, Central and South American disjunct pattern is not unusual in angiosperms, and in fact has been investigated in several studies (e.g., Wen et al., 2010), the species-rich family Asteraceae provides a good system to probe the biogeographic pattern of taxa with distinct distribution patterns around the globe (Panero and Funk, 2008). Specifically members of Adenocaulon may provide an excellent opportunity to examine species adaptation in a recent diversification. Fossil and phylogenetic evidence demonstrates that Asteraceae may be of South American origin and underwent a major diversification event, followed by an African explosion (DeVore and Stuessy, 1995; Funk et al., 2005; Barreda et al., 2010; Stuessy, 2010) and subsequently moved to various continents (Nie et al., 2016). However, this model does not apply to the small genus Adenocaulon.

The goal of this study is to elucidate the phylogenetic relationships within the genus and to infer its historical biogeography utilizing chloroplast and nuclear DNA in conjunction with fossil evidence. Knowing details about how and when the genus originated and dispersed can provide insight into understanding the plausible biogeographic affinity between East Asia and America. Our study would shed new lights on the formation of intercontinental disjunction of plants.

## MATERIALS AND METHODS

## Taxa Sampling and Data Collection

Ten samples (accessions) were used for each of the 5 species of Adenocaulon (**Table 1**). Based on the results from previous phylogenetic analyses of Asteraceae (Panero and Funk, 2008; Ortiz et al., 2009), additional samples were used from eight genera that are closely related to Adenocaulon (Gerbera, Leibnitzia, Chaptalia, Trichocline, Brachyclados, Pachylaena, Mutisia, and Chaetanthera of the tribe Mutisieae of the subfamily Mutisioideae) and 15 distantly related genera of the subfamily (Trixis, Dolichlasium, Jungia, Perezia, Nassauvia, Acourtia, Leucheria, Proustia, and Lophopappus in the tribe Nassauvieae of the subfamily Mutisioideae; Onoseris, Urmenetea, Lycoseris, Aphyllocladus, Gypothamnium, and Plazia from the tribe


*(Continued)*

TABLE 1 | Species names, localities, voucher deposition,

 and GenBank accession numbers used in the present study.


Onoserideae of the subfamily Mutisioideae). We also included nine species from Wunderlichioideae and five species from Barnadesioideae. We used Acicarpha spathulata (Calyceraceae) to serve as the outgroup, following Panero and Funk (2008) (**Table 1**).

## DNA Extraction, Amplification, and Sequencing

Total genomic DNA was isolated from silica-dried leaf materials and herbarium specimens using a Universal Genomic DNA Extraction Kit (Takara, Dalian, China). Nine gene regions, including eight plastid DNA regions (ndhF, matK, rpoB, rbcL, trnL-F, rpl32-trnL, rpl32-ndhF, and psbA-trnH) and the nuclear ITS, were employed. According to Shaw et al. (2007) and Timme et al. (2007), both rpl32-trnL and rpl32-ndhF have faster molecular evolution rate than other chloroplast regions and can provide adequate resolution in phylogenetic reconstructions of Asteraceae. The ITS data were generated from single amplifications using primers ITS4 and ITS5 (White et al., 1990). Primers used for amplification and sequencing were tabe and tabf for trnL-F region (Taberlet et al., 1991), Z1 and 1204R for rbcL (Zurawski et al., 1981), and psbA\_F and trnH\_R for psbA-trnH intergenic spacer (Sang et al., 1997). The primer sets for ndhF, matK, and rpoB regions were from Panero and Funk (2008), while the primers for rpl32-trnL and rpl32-ndhF were from Shaw et al. (2007). Amplified products were purified using a Qiaquick gel extraction kit (Qiagen, Inc., Valencia, California, USA) and sequenced in both directions by an ABI 3730 automated sequencer (Applied Biosystems, Foster City, California, USA). The resulting sequences were edited using Sequencher (version 4.1.4) and aligned with MUSCLE (version 3.6; Edgar, 2004), followed by manual adjustment in Se-Al (v2.0a11; Rambaut, 2002). All sequences generated for this study were deposited in GenBank under accession numbers shown in **Table 1**.

## Phylogenetic Analysis

The nuclear and the plastid datasets were analyzed both separately and simultaneously using maximum parsimony (MP), Bayesian inference (BI), and maximum likelihood (ML). The topologies of the data sets were compared with each other to detect any incongruence. Because no significant incongruence was observed, we chose to combine the nrDNA and cpDNA data sets. For the maximum parsimony analysis we used PAUP/(version 4.0b10; Swofford, 2003). All characters were weighted equally and unordered. Heuristic searches were conducted using 100 random-taxon-addition replicates with tree bisection-reconnection (TBR) branch swapping, using MulTrees option in effect, and a maximum of 10,000 trees. Bootstrap analyses (1,000 pseudoreplicates) were conducted to examine the relative level of support for individual clades on the cladograms of each search (Felsenstein, 1985).

Models of nucleotide substitution were selected based on the Akaike Information Criterion (AIC) as determined by MrModelTest (version 2.3; Nylander, 2004). Maximumlikelihood searches and bootstrap analyses were performed on the XSEDE online computing cluster accessed via the CIPRES Science Gateway (Miller et al., 2010) using RAxML-HPC2

Bayesian inference was conducted using MrBayes (version 3.2.2, Ronquist and Huelsenbeck, 2003). The Markov chain Monte Carlo (MCMC) algorithm was run for 10,000,000 generations with one cold and three heated chains, starting from random trees and sampling one out of every 1,000 generations. The burn-in and convergence diagnostics were graphically assessed using AWTY (Nylander et al., 2008). The burn-in trees were excluded, and the remaining trees were assumed to be representative of the posterior probability (PP) distribution.

## Divergence Time Estimation

The combined nuclear and plastid sequences were used to estimate the divergence time of Adenocaulon using BEAST (version 1.8.0; Drummond and Rambaut, 2007). BEAST employs a Bayesian MCMC approach to co-estimate topology, substitution rates and node ages (Drummond et al., 2002). BEAUti was used to set criteria for the analysis, in which we applied a general time reversible (GTR) nucleotide-substitution model with Gamma + Invariant sites, gamma shape distribution (with four categories) and proportion of invariant sites. A Yule tree prior model was implemented in the analysis, with rate variation across branches assumed to be uncorrelated and lognormally distributed (Drummond et al., 2006). Posterior distributions of parameters were approximated using two independent MCMC analyses of 50,000,000 generations (sampling once for every 5,000 generations) after a 10% burn-in for each. Convergence of the chains was checked using Tracer (version 1.5; Rambaut and Drummond, 2007).

For the fossil calibrations, we followed the strategy adopted by Nie et al. (2016). The first calibration point was the Raiguenrayun cura Barreda, Katinas, Passalia & Palazzesi capitulescence from the Patagonia described as belonging to the crown Asteraceae and dated by radiometric methods at 47.5 Ma (von Nickisch-Rosenegk et al., 1999). This fossil was used as a minimum age constraint for the split between Barnadesioideae and the rest of family and with an applied lognormal prior distribution with an offset of 47.5 Ma, a mean of 1.0 and a standard deviation of 0.7. As a minimum age constraint for the crown group of the subfamily Barnadesioideae, we used the fossil pollen Quilembaypollis sp., from the latest Oligocene/earliest Miocene, dated at 23 Ma (Chiang et al., 1998). Following Bergh and Linder (2009), the 95% confidence interval for this prior lies between 16.9 and 44.1 Ma with the mean at 22.3 Ma.

## Biogeographic Analyses

We compiled distribution data for the Adenocaulon species and assigned the included taxa to respective ranges. Considering the areas of endemism in Adenocaulon and tectonic history of continents, we used four areas for the Adenocaulon taxa: (A) South America, (B) Central America, (C) North America, and (D) East Asia. Each sample was assigned to its respective area according to its contemporary distribution range.

Biogeographic inferences were obtained by applying both statistical dispersal–vicariance analysis (S-DIVA) and dispersalextinction-cladogenesis (DEC) model calculated by Lagrange (Ree and Smith, 2008) implemented in RASP (version 3.0) using the default settings (Yu et al., 2010). A subset of 1,000 randomly selected trees from the posterior distribution output of BEAST was used and the maximum number of individual unit areas was set to two. The probability of dispersal between areas was maintained as equal.

## RESULTS

## Phylogenetic Analyses

The incongruence length difference and Shimodaira and Hasegawa (SH) tests failed to reveal significant incongruence, allowing the individual datasets to be combined. There were 12,416 characters of combined data set, of which 1,571 were variable, and 1,266 parsimony informative sites. In general, the combined data matrix provided greater resolution and stronger support for phylogenetic relationships than did the individual datasets (**Figure 2**).

The results from MP and BI analysis, as well as ML analysis, strongly support Adenocaulon being a monophyletic clade (PP = 100, BP = 100, PL = 100). Within Adenocaulon, the Central American species, A. lyratum, diverged first, followed by A. chilense being sister to a clade including the remaining species from North America and East Asia albeit with low support (PP = 64, BP = 54, PL = 85). The two species found in East Asia, A. himalaicum and A. nepalense, form a clade with high support values (PP = 100, BP = 93, PL =98), which is sister to A. bicolor (North America) (**Figure 2**).

## Biogeographic Analysis

The chronogram and results of divergence-time estimation based on the nrDNA and cpDNA data showed that the crown age of Adenocaulon was estimated at 9.05 Ma with a 95% highest posterior density (HPD) of 4.73–14.99 Ma. The stem age, that is the divergence time between Adenocaulon and its close relatives, was estimated at 14.74 Ma (with 95% HPD: 9.6–23.7 Ma). The split between the South American species A. chilense from the remaining species of Adenocaulon was at 6.92 Ma (95%HPD: 3.55–12.05 Ma), whilst the disjunction between North American A. bicolor and East Asian species at 4.79 Ma (95% HPD: 2.11– 8.68) (**Figure 3**). S-DIVA analysis demonstrated that Adenocaulon may have originated in the tropical America and then dispersed to North America and East Asia (**Figure 3**).

## DISCUSSION

## Phylogenetic Relationship of Adenocaulon with Other Genera within Mutisieae

Our study showed that Adenocaulon is monophyletic within the Mutisieae (PP = 100, BP = 100, PL = 100) and is sister to a clade including Gerbera, Leibnitzia, Chaptalia, Trichocline, and Brachyclados (**Figure 2**). The genus is characterized by

a distinct set of synapomorphy including basal constriction of anther appendages, small anthers, disciform capitula, and typical glandular achenes without pappus (Bittmann, 1990).

The placement of Adenocaulon has long been disputed among numerous taxonomists, going as far as calling it an anomalous genus, and being unplaced in the tribe (Bremer,

nodes are bootstrap values obtained from MP analysis and posterior probabilities support values from ML analysis. "\*" indicates full support for the respective analysis.

1994). Based on our broad analysis of the early diverging clades of Asteraceae, Adenocaulon is confirmed to belong to the tribe Mutisieae (Mutisioideae) and to be sister to the Gerbera-complex (PP = 100, BP = 79, PL = 76) (**Figure 2**).

Our results indicate that the Central American A. lyratum diverged early. This species differs from the other four species in a stem winged by decurrently leaf bases, lyrate-pinnatifid sessile leaves and distinct club-shaped achene. Representatives of A. chilense from southern Chile and Argentina show petiolate entire leaves in a rosulate order. The pistillate flowers have a distinct bilabiate corolla and the scabrate pollen exhibit two layers of columellae with a larger inner one. The other three species from the Northern Hemisphere form a clade with high support. They share common features, e.g., leaves concentrated at the base, leaves triangular with cordate leaf base, heart- or kidney-shaped, with dentate or sinuate margin. Within this clade, the North American A. bicolor is sister to a clade with the two Asian species A. bicolor and A. himalaicum; these two species can be easily distinguished from each other by the glandular hairs: the former being pin-shaped, occurring on both vegetative and synflorescence branches, whereas the latter is nail-shaped but restricted to the synflorescence branches (Bittmann, 1990). A. nepalense is a recently described new species (Bittmann, 1990) endemic to Bhutan, Nepal, and India, and having distinct features such as winged by decurrently leaf-based stems and the bearing whip hairs achenes.

## The Origin and Early Dispersal of Adenocaulon

Our biogeographic reconstruction suggests that Adenocaulon may have originated in Central America and then dispersed to North America and South America during the Miocene (14.74 Ma, 95% HPD: 9.6–23.7) (**Figure 3**). The dispersal of biota between Central America, North America and South America (commonly called the Great American Biotic Interchange) is one of the most significant events in neotropical biogeography. The Great American Biotic Interchange was caused by the closure of the Isthmus of Panama ca. 3 Ma (Burnham and Graham, 1999; Cody et al., 2010), resulting in the exchange of plants and animals between the two once separated continents. However, some authors suggest that an earlier land connection made possible the exchange of biota within the Americas. A land bridge may have existed between North and South Americas 3–7 Ma (Bermingham and Martin, 1998). Moreover, some plants could have successfully crossed the Isthmus of Panama before it closed (Cody et al., 2010; Bacon et al., 2015). Here, we propose that Adenocaulon may have dispersed into South America from Central America during the Miocene via the land bridge link. Hypotheses of long distance dispersal events have gained popularity in the last decade as a way of explaining many disjunct distributions of plant lineages (Givnish et al., 2004; Popp et al., 2011; Nie et al., 2013).

Dispersal of Adenocaulon may have been enhanced by birds due to its special fruit structures. The prominent glandular hairs on the cypselae are very sticky and cling readily to fabrics, fur, and feathers. The intercontinental disjunction within the range of Adenocaulon is possibly a result of dispersal of cypselae on the feathers of birds. The bird dispersal pattern may also account for the intracontinental disjunction of A. bicolor (in the Pacific Northwest, the Black Hills, and the northern Great Lakes region). Thus, we do not postulate vicariance as a factor that shaped the disjunct pattern of Adenocaulon.

## The Intercontinental Dispersal within Northern Hemisphere

The intercontinental dispersal of plants in the Northern Hemisphere has been a focus of biogeography for a long time (Donoghue et al., 2001; Donoghue and Smith, 2004). During the Cenozoic, two routes for plant dispersal between the New World and the Old World were the Bering Land Bridge and the North Atlantic Land Bridge (Milne and Abbott, 2002; Wen et al., 2010). These two links have played significant roles in the formation of the modern flora of the Northern Hemisphere but happened in different geologic times (Tiffney, 1985a,b). The Bering Land Bridge is assumed to have existed from the Paleocene to Miocene (Tiffney and Manchester, 2001) but later this dispersal route was blocked due to global cooling (Milne, 2006). We suggest that Adenocaulon most likely migrated from North America to East Asia via the Bering Land Bridge. The estimated divergence time between the North American and East Asian groups is ca. 11.91 Ma, which coincided the cooling of Northern Hemisphere in which the continuous distribution of the Boreotropical flora was disrupted (Tiffney and Manchester, 2001). Our results suggest Adenocaulon may have originated in Central America, dispersed into North America and further into East Asia via the Bering Land Bridge, and finally reached the Himalayas.

## Role of Fruits in Adenocaulon Dispersibility

Our results suggest that the biogeographic history of Adenocaulon species have been influenced by dispersal events, which could be facilitated by their unusual floral features. The pappus and the fruit trichome can prevent desiccation, and other characteristics of the family lack in Adenocaulon. Instead, the fruits of Adenocaulon are covered by glands with a sticky secretion. An umbrella-like "infructescence" emerges from the otherwise short plants usually embedded in other under-story vascular plants or mosses, and exposes the sticky fruits. Pervious researches indicated that viscid fruits can adhere by bird feathers (Carlquist, 1967, 1983). Also, the possibility that birds could have been the dispersal vectors for Adenocaulon fruits was proposed by Bittmann (1990) to explain some of the extant distributions of the genus. Alternatively, successful over-water dispersal events could be postulated for Adenocaulon; the oily fruit cover is suggested to aid in flotation, as was the case for other families. Thus, the absence of pappus in Adenocaulon does not mean that this genus lacks mechanisms for dispersal. The sticky fruits, with the "infructescence," can act as an alternative dispersal mechanism involving several potential dispersal vectors.

## AUTHOR CONTRIBUTIONS

HS: conceived this research project; TD, YC, HW, XZ, SV, ZY, and HQ: collected the data; TD and YC: analyzed the data; TD, YC, SV, HW, and HS: led the writing.

## ACKNOWLEDGMENTS

This study was supported by the Major Program of National Natural Science Foundation of China (31590823) and the National Key R&D Program of China (2017YFC0505200) to HS, and the National Natural Science Foundation of China (31700165) and CAS "Light of West China" Program to TD. We thank Dr. Jacob Landis of University of California, Riverside in English editing and Dr. Rinchen Yangzom of National Herbarium, National Biodiversity Center of Bhutan for providing leaf materials.

## REFERENCES


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Deng, Chen, Wang, Zhang, Volis, Yusupov, Qian and Sun. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Climatic Controls on C<sup>4</sup> Grassland Distributions During the Neogene: A Model-Data Comparison

David L. Fox <sup>1</sup> \*, Stephanie Pau<sup>2</sup> , Lyla Taylor <sup>3</sup> , Caroline A. E. Strömberg<sup>4</sup> , Colin P. Osborne<sup>3</sup> , Catherine Bradshaw5†, Stephen Conn5†, David J. Beerling<sup>3</sup> and Christopher J. Still <sup>6</sup>

#### Edited by:

Oana Moldovan, Emil Racovita Institute of Speleology, Romania

#### Reviewed by:

David Nelson, University of Maryland Center for Environmental Science (UMCES), United States Lydie M. Dupont, Zentrum für Marine Umweltwissenschaften, Universität Bremen, Germany

> \*Correspondence: David L. Fox dlfox@umn.edu

#### †Present Address:

Catherine Bradshaw, Met Office Hadley Centre, Exeter, United Kingdom Stephen Conn, School of Earth and Ocean Sciences, Cardiff University, Cardiff, United Kingdom

#### Specialty section:

This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution

Received: 21 May 2018 Accepted: 04 September 2018 Published: 24 September 2018

#### Citation:

Fox DL, Pau S, Taylor L, Strömberg CAE, Osborne CP, Bradshaw C, Conn S, Beerling DJ and Still CJ (2018) Climatic Controls on C4 Grassland Distributions During the Neogene: A Model-Data Comparison. Front. Ecol. Evol. 6:147. doi: 10.3389/fevo.2018.00147 <sup>1</sup> Department of Earth Sciences, University of Minnesota, Minneapolis, MN, United States, <sup>2</sup> Department of Geography, Florida State University, Tallahassee, FL, United States, <sup>3</sup> Department of Animal and Plant Sciences, University of Sheffield, Sheffield, United Kingdom, <sup>4</sup> Department of Biology, University of Washington, Seattle, WA, United States, <sup>5</sup> Bristol Research Initiative for the Dynamic Global Environment, School of Geographical Sciences, University of Bristol, Bristol, United Kingdom, <sup>6</sup> Department of Forest Ecosystems and Society, Oregon State University, Corvallis, OR, United States

Grasslands dominated by taxa using the C<sup>4</sup> photosynthetic pathway first developed on several continents during the Neogene and Quaternary, long after C<sup>4</sup> photosynthesis first evolved among grasses. The histories of these ecosystems are relatively well-documented in the geological record from stable carbon isotope measurements (of fossil vertebrate herbivores and paleosols) and the plant microfossil record (pollen and/or phytolith assemblages). The distinct biogeography and ecophysiology of modern C<sup>3</sup> and C<sup>4</sup> grasses have led to hypotheses explaining the origins of C<sup>4</sup> grasslands in terms of long-term changes in the Earth system, such as increased aridity and decreasing atmospheric pCO2. However, quantitative proxies for key abiotic drivers of these hypotheses (e.g., temperature, precipitation, pCO2) are still in development, not yet widely applied at the continental or global scale or throughout the late Cenozoic, and/or remain contentious. Testing these hypotheses globally therefore remains difficult. To understand better the potential links between changes in the Earth system and the origin of C<sup>4</sup> grasslands, we undertook a global-scale comparison between observational records of C<sup>4</sup> plant abundances in Miocene and Pliocene localities compiled from the literature and three increasingly complex models of C<sup>4</sup> physiology, dominance, and abundance. The literature compilation comprises >2,600 δ <sup>13</sup>C-values each of fossil terrestrial vertebrates and of paleosol carbonates, which we interpret as primarily proxies for the abundance of C<sup>4</sup> grasses, based on the modern contribution of C<sup>4</sup> grasses to terrestrial net primary productivity. We forced the vegetation models with simulated monthly climates from the HadCM3 family of coupled ocean-atmosphere general circulation models (OAGCMs) over a range of pCO2-values for each epoch to model C<sup>4</sup> dominance or abundance in grid cells as: (1) months per year exceeding the temperature at which net carbon assimilation is greater for C<sup>4</sup> than C<sup>3</sup> photosynthesis (crossover temperature model); (2) the number of months per year exceeding the crossover temperature and having sufficient precipitation for growth (≥25 mm/month; Collatz model); and (3) the Sheffield Dynamic Global Vegetation Model (SDGVM), which models multiple plant functional types (PFTs) (C<sup>3</sup> and C<sup>4</sup> grasses, evergreen, and deciduous trees). Model-data agreement is generally weak, although statistically significant for many

**44**

comparisons, suggesting that regional to local ecological interactions, continent-specific plant evolutionary histories, and/or regional to local climatic conditions not represented in global scale OAGCMs may have been equally strong or stronger in driving the evolution of C<sup>4</sup> grasslands as global changes in the Earth system such as decreases in atmospheric pCO<sup>2</sup> and late Cenozoic global cooling and/or aridification.

Keywords: Miocene, Pliocene, C4 grasses, carbon isotopes, model-data comparison, vegetation models

## INTRODUCTION

Plants using the C<sup>4</sup> photosynthetic pathway comprise only a small fraction of the species diversity of the modern global flora, and about half of these extant C<sup>4</sup> species are grasses. Despite the limited number of C<sup>4</sup> species, tropical, and subtropical ecosystems dominated by C<sup>4</sup> grasses account for almost one-quarter of modern terrestrial net primary productivity (Still et al., 2003; Sage, 2004; Sage et al., 2011). Molecular phylogenies indicate that the C<sup>4</sup> syndrome evolved multiple times independently in grasses (Kellogg, 1999; GPWGII—Grass Phylogeny Working Group II, 2012), and the calibration of these phylogenies to time using the fossil record suggests that the earliest C<sup>4</sup> grasses had evolved by the Oligocene epoch (33– 25 Ma), with successive independent origins of C<sup>4</sup> anatomy and biochemistry continuing throughout the Neogene Period (25– 2.6 Ma; Christin et al., 2008, 2014; Spriggs et al., 2014).

The increasingly comprehensive species sampling of grass phylogenies has dramatically improved our ability to generate testable hypotheses about the timing of trait evolution and the diversification of C<sup>4</sup> lineages, but cannot reveal the temporal pattern of increasing biomass and the ecological importance of C<sup>4</sup> grasses in local and regional ecosystems. Many grasslands are dominated locally by only a few species in terms of frequency, abundance, and biomass (Smith and Knapp, 2003; Edwards et al., 2010; Griffith et al., 2015), such that species diversification over time does not necessarily predict the ecological success of C<sup>4</sup> grasses generally. However, the geological record provides evidence about both when grasslands emerged and when they became dominated by C<sup>4</sup> species.

The plant fossil record, particularly phytoliths (microscopic intra- and extracellular silica bodies diagnostic to specific clades of grasses), indicates that grass-dominated habitats first appeared on some continents by the late Oligocene to early Miocene epochs. Complementary evidence based on carbon isotope records from fossil teeth of mammalian herbivores and calcareous paleosols (fossilized soils that preserve soil-derived carbonate) shows that C4-dominated grasslands first appeared on various continents during the late Miocene to early Pleistocene (Jacobs et al., 1999; Tipple and Pagani, 2007; Edwards et al., 2010; Strömberg, 2011; Fox et al., 2012a). However, various terrestrial isotopic records (e.g., Kleinert and Strecker, 2001), as well as marine records of terrigenous materials (e.g., Feakins et al., 2013; Hoetzel et al., 2013), indicate that the temporal pattern of the emergence of C4-dominated grasslands was not a simple, monotonic increase in C<sup>4</sup> biomass in some regions. Instead, some records (e.g., the terrigenous signal in marine sediments of northeast African vegetation during the Neogene; Feakins et al., 2013) indicate fluctuations through time in the C3:C<sup>4</sup> ratio at regional scales prior to the emergence of modern, C4-dominated grasslands. Thus, the ecological dominance of C<sup>4</sup> grasses in many grassland biomes today post-dates the origins of most C<sup>4</sup> grass clades by many millions of years and reflects a complex history of ecosystem assembly. Despite the attention of a diverse community of systematists, ecologists, paleontologists, and geochemists over many years, the specific combinations of environmental and climatic changes responsible for the ecological success of C<sup>4</sup> grasses over the Neogene and Quaternary both globally and in specific regions have remained elusive and debated (see discussions in, e.g., Edwards et al., 2010; Scheiter et al., 2012; Bond, 2015).

C<sup>4</sup> photosynthesis describes a suite of anatomical and biochemical characteristics that together act as a CO<sup>2</sup> pump to reduce photorespiration, thereby increasing the quantum yield, a measure of how efficiently plants fix CO<sup>2</sup> relative to the amount of absorbed photosynthetically active radiation (APAR) received during light-limited photosynthesis (Ehleringer and Björkman, 1977; Collatz et al., 1998). Photorespiration rates are higher at warmer temperaturesfor a given atmospheric pCO2, and at lower pCO<sup>2</sup> for a given temperature (Ehleringer et al., 1997; Sage, 2004). The dominance of modern C<sup>4</sup> grasses in warm to hot and mesic to arid regions (e.g., Teeri and Stowe, 1976; Hattersley, 1983; von Fischer et al., 2008) has led to a long-held assumption that C<sup>4</sup> photosynthesis evolved as an adaptation to warm and dry climates. However, ancestral state reconstructions of the climatic origins of C<sup>4</sup> lineages have called this assumption into question, by showing that most C<sup>4</sup> grass lineages derive from C<sup>3</sup> ancestors that already inhabited warm climates (Edwards and Still, 2008; Edwards and Smith, 2010). From these origins in warm climates, C<sup>4</sup> photosynthesis enabled diversification into cooler and drier climates (Osborne and Freckleton, 2008; Watcharamongkol et al., 2018).

Some paleoclimate records associated with estimates of C<sup>4</sup> biomass, such as compound-specific isotope analyses of terrestrial leaf waxes in Pliocene marine sediments off of West Africa (Kuechler et al., 2018), are consistent with a relationship between aridity and C<sup>4</sup> biomass during a time with atmospheric pCO<sup>2</sup> ≥300 ppmV (i.e., above pre-industrial values; Foster et al., 2017). However, similar records from Pleistocene marine sediments off of North Africa (e.g., Kuechler et al., 2013) indicate an inverse relationship between regional aridity in northern Africa and C<sup>4</sup> abundance against a background of low, but fluctuating glacial-interglacial pCO<sup>2</sup> (see also Urban et al., 2015). These contrasting records suggest that the role of aridity in controlling the abundance of C<sup>4</sup> plants in regional ecosystems may be strongly modulated by factors other than CO2. These observations are broadly consistent with our understanding of C<sup>4</sup> grass biogeography today. For example, independently evolved C<sup>4</sup> grass lineages show contrasting relationships with rainfall (Visser et al., 2012, 2014). In addition, very arid areas are rarely dominated by C<sup>4</sup> species, whereas particular C<sup>4</sup> grass clades may dominate humid regions if the dry season is sufficient to allow frequent fires (Bond, 2015).

Decreasing atmospheric pCO<sup>2</sup> during the Cenozoic was initially proposed as an explanation for the appearance of C4-dominated grasslands beginning in the late Miocene (e.g., Cerling et al., 1997; Ehleringer et al., 1997). Subsequent work challenged this assertion, demonstrating that atmospheric pCO<sup>2</sup> first decreased to levels that could favor C<sup>4</sup> over C<sup>3</sup> photosynthesis given prevailing climatic conditions during the Oligocene (Pagani et al., 1999; Beerling and Royer, 2011; Foster et al., 2017), not the late Miocene (although note that some CO2 proxy records indicates a rise and subsequent fall in pCO<sup>2</sup> during the middle-late Miocene; e.g., Kürschner et al., 2008; Bolton and Stoll, 2013). The broad coincidence in the timing of declining atmospheric pCO<sup>2</sup> and the proposed, earliest evolution of C<sup>4</sup> photosynthesis in grasses has therefore pointed to a role for low pCO<sup>2</sup> in the evolution of C<sup>4</sup> photosynthesis (e.g., Christin et al., 2008; Spriggs et al., 2014). In contrast, the asynchronous appearance of C4-dominated grasslands in different continents from the late Miocene to the early Pleistocene (Fox and Koch, 2003; Edwards et al., 2010) suggests that additional environmental factors, in combinations specific to each region, were of equal or larger importance for increasing the dominance of C<sup>4</sup> grasses in many ecosystems. These factors could include both biotic and abiotic influences such as disturbance by herbivores, changes in fire frequency and intensity, differential patterns of recovery by C<sup>3</sup> trees and grasses and C<sup>4</sup> grasses linked to traits other than photosynthetic pathway (i.e., evolutionary history of different lineages), soil type and characteristics, and region-specific patterns of climate change during the Neogene and Quaternary (e.g., Osborne, 2008; Edwards et al., 2010; Bond, 2015; Griffith et al., 2015; Charles-Dominique et al., 2016). The importance of declining CO<sup>2</sup> as a single factor in promoting C<sup>4</sup> dominance has also been called into question by the results of a 20 year free-air CO<sup>2</sup> enrichment experiment that showed an eventual enhancement in C<sup>4</sup> biomass relative to C<sup>3</sup> biomass growing under experimentally increased CO<sup>2</sup> (Reich et al., 2018). Water cycle feedbacks induced by enhanced CO<sup>2</sup> have also been shown experimentally to cause C<sup>4</sup> grasses to outperform C<sup>3</sup> grasses, with greater water-use efficiency by C<sup>4</sup> grasses driven by interacting effects of higher CO<sup>2</sup> and warming on both soil moisture and plant photosynthetic responses as the proposed mechanism (Morgan et al., 2011).

Fire in particular is recognized as a critical factor in the abundance of C<sup>4</sup> grasses in many modern ecosystems (Bond, 2015), and fire has been proposed as playing an important role in the emergence of C4-dominated grasslands during the late Neogene (Keeley and Rundel, 2005; Scheiter et al., 2012). The importance of fire for C4-dominated grasslands has support from micro-charcoal counts in both Pliocene (Hoetzel et al., 2013) and Holocene (Dupont and Schefuss, 2018) marine sediments from cores off of West Africa, but the terrestrial record of charcoal is still too incomplete (e.g., Scott, 2000) to allow for direct comparisons of charcoal occurrences with either paleobotanical or vertebrate and paleosol carbon isotope proxies documenting the emergence of C4-dominated ecosystems.

A wide range of continental paleoclimate proxies derived from paleontological and geological materials can provide critical information on the abiotic conditions associated with the increase of C<sup>4</sup> biomass in local and regional ecosystems. Some of these, such as the oxygen isotope composition of fossil teeth of mammalian herbivores (e.g., Passey et al., 2002) or paleosol carbonates (e.g., Fox et al., 2012b) are controlled by both temperature and precipitation, so cannot generally be deconstructed to provide unambiguous quantitative estimates of paleoclimatic conditions. Other proxies, such as the relative tooth crown heights of mammalian herbivores (Eronen et al., 2010) or weathering indices based on elemental ratios of paleosols (Sheldon and Tabor, 2009; Stinchcomb et al., 2016), can provide quantitative estimates of mean annual temperature and precipitation. However, such proxies, as well as newer methods such as carbonate clumped isotope paleothermometry (Ghosh et al., 2006; Eiler, 2011), either rely on assumptions about faunal trait evolution that have been called into question (e.g., MacFadden et al., 1999; Strömberg, 2006), or have not yet been applied widely enough to allow for compilations at the continental or global scale with sufficient coverage to allow for meaningful analyses of empirical paleoclimate data in relation to records of the Neogene-to-Quaternary rise in global C<sup>4</sup> biomass. In contrast, global scale, ocean-atmosphere general circulation models (OAGCMs) coupled with vegetation modeling provide a means to assess the record of C<sup>4</sup> biomass on all continents in relation to possible past climatic conditions that ultimately can be tested with paleoclimate proxy data (e.g., see Lunt et al., 2007 for a comparison of late Oligocene and pre-industrial boundary conditions).

In this paper, we undertake a global-scale comparison between a literature compilation of δ <sup>13</sup>C values of mineralized tissues of fossil herbivorous vertebrates (mammal teeth, ratite eggshell) and of paleosol carbonates from Miocene and Pliocene localities as proxies for the abundance of C<sup>4</sup> grasses, and three increasingly complex models of C<sup>4</sup> dominance and abundance driven by output from a global OAGCM. Specifically, we use climate output from the HadCM3LBL-M2.1 and HadCM3BL-M1 OAGCMs (Valdes et al., 2017) over a range of pCO<sup>2</sup> values for the Miocene (Bradshaw et al., 2012) and Pliocene (Bragg et al., 2012; Conn, 2012) to model C<sup>4</sup> dominance or abundance in grid cells as: (1) months per year exceeding the temperature at which net assimilation is greater for C<sup>4</sup> than C<sup>3</sup> photosynthesis (crossover temperature model; Ehleringer et al., 1997); (2) the number of months per year exceeding the crossover temperature and having sufficient precipitation for grass growth (≥25 mm/month in months that meets the crossover temperature criterion; Collatz model; Collatz et al., 1998), and (3) predictions of various measures of C<sup>4</sup> biomass from the Sheffield Dynamic Global Vegetation Model (SDGVM; Woodward and Lomas, 2004), output from which includes a range of vegetation parameters of multiple plant functional types (PFTs), including C<sup>3</sup> and C<sup>4</sup> grasses. To account for the uncertainty in estimates and temporal variation in CO<sup>2</sup> levels suggested by proxy data for this time interval (see Beerling and Royer, 2011), we used a range of pCO<sup>2</sup> values in each analysis. Statistical comparisons of the isotopic databases with the paleoclimate and vegetation model outputs allow us to assess data-model agreement. We propose that close similarity between proxy data and a particular model will provide support for the explanatory factors inherent in the model playing an important role in C<sup>4</sup> dominance globally (or regionally). A mismatch, on the other hand, will suggest that the factors accounted for in the model are not sufficient for explaining the patterns seen in the data, although alternative reasons are possible. Specifically, we discuss the limitations and opportunities with current data and model resolution for inferring the role of abiotic factors in the evolution of the modern C<sup>4</sup> grasslands over the Neogene and Quaternary.

## CARBON ISOTOPE PROXIES FOR C<sup>4</sup> BIOMASS

The basis for using the stable carbon isotope composition (δ <sup>13</sup>C) of carbonate-bearing minerals of vertebrate herbivores and authigenic soil carbonates as proxies for C<sup>4</sup> biomass is the difference in isotope discrimination by C<sup>3</sup> and C<sup>4</sup> plants during fixation of CO2. Due to physiological differences, C<sup>3</sup> plants discriminate more strongly against <sup>13</sup>CO<sup>2</sup> than do C<sup>4</sup> plants (O'Leary, 1981), and, consequently, tissues of C<sup>3</sup> plants have much lower (more negative) δ <sup>13</sup>C values than do C<sup>4</sup> plants. Here, we treat the δ <sup>13</sup>C values of modern plants reported by Passey et al. (2002) as end-members, whereby modern C<sup>3</sup> plants have mean δ <sup>13</sup>C (±1 s.d.) of −27.4 ± 1.6% and modern C<sup>4</sup> monocots have mean δ <sup>13</sup>C of −12.7 ± 1.1%. Passey et al. (2002) assumed an average δ <sup>13</sup>C value for atmospheric CO<sup>2</sup> for these plant data of −8.0%, which allowed estimation of apparent fractionation factors during fixation of CO<sup>2</sup> by C<sup>3</sup> and C<sup>4</sup> photosynthesis of −19.6 and −4.7%, respectively.

Authigenic carbonate precipitates in soils in isotopic equilibrium with the CO<sup>2</sup> that is dissolved in soil water (Cerling, 1984). For soils with moderate to high respiration rates and at modern atmospheric pCO2, all CO<sup>2</sup> below about 30 cm in the profile is derived from the standing plant biomass or respiration of plant-derived components in the soil, and has the carbon isotopic signature of the standing biomass (Cerling et al., 1991). Relative to atmospheric CO2, pedogenic carbonate is enriched in <sup>13</sup>C by ca. 14–17% due to faster diffusion of <sup>12</sup>CO<sup>2</sup> from the soil to the atmosphere and the temperature-dependent fractionation of carbonate precipitation over the range of typical soil temperatures (Cerling et al., 1991). Paleosol carbonates form on timescales of 10<sup>2</sup> -10<sup>5</sup> years, depending on sedimentation rates and landscape position (Gile et al., 1966; Birkeland, 1999); however, they sample the landscape on very localized spatial scales assuming minimal horizontal advection of CO<sup>2</sup> in the subsurface (i.e., by soil water flow). Importantly, vegetation, soil type and properties, and presence or preservation of carbonate can all vary on small spatial scales with landscape position (Zamanian et al., 2016). Thus, the δ <sup>13</sup>C of paleosol carbonates preserves a time-integrated signature of plant biomass that can be used to estimate the C3:C<sup>4</sup> ratio of the local ecosystem on relatively long timescales but small spatial scales.

The δ <sup>13</sup>C values of vertebrate herbivore tissues record the δ <sup>13</sup>C of diet with tissue-specific fractionation factors or apparent enrichment factors (in the case of tissues such as enamel for which the biosynthetic reactions are not reversible, DeNiro and Epstein, 1978; Cerling and Harris, 1999; Passey et al., 2005). For structurally bound carbonate in the hydroxyapatite of tooth enamel of large-bodied mammalian herbivores, the apparent enrichment factor relative to diet is +14.1 ± 0.5% (Cerling and Harris, 1999). The isotopic composition of mammalian tooth enamel is generally considered to be resistant to alteration during fossilization (Wang and Cerling, 1994; Koch et al., 1997; Zazzo et al., 2004), so that fossil teeth of herbivorous mammals faithfully record the dietary proportions of C<sup>3</sup> and C<sup>4</sup> plants during the interval of tooth formation (Kohn and Cerling, 2002). Additionally, Griffith et al. (2017) showed that the δ <sup>13</sup>C-values of grazer tissue (bison, mammoth) faithfully records the C3:C<sup>4</sup> ratio of grasslands across larger spatial scales. The carbonate of ratite eggshell also records the δ <sup>13</sup>C of the diet of the hen while the egg is mineralizing (Schirnding et al., 1982), and eggshell δ <sup>13</sup>C values also appear to be resistant to alteration during fossilization (Stern et al., 1994; Miller et al., 2005). Fossil deposits are subject to time averaging, such that data for fossils from individual localities can represent a wide range of ages (Behrensmeyer et al., 2000). However, individual teeth and eggshells only represent a few months to years of mineralization; thus, the temporal resolution of teeth is higher than that for paleosols. On the other hand, the animals did not necessarily live where the fossils accumulated (Behrensmeyer et al., 2000), and mammalian home range generally scales positively with body size (Harestad and Bunnel, 1979). For example, a largebodied mammal like a moose has a home range of ca. 1,500 ha, but a smaller bodied species such as a peccary has a home range of only ca. 135 ha (data from Harestad and Bunnel, 1979). Teeth of large mammalian herbivores, which have been most commonly used for paleovegetation reconstruction, therefore integrate the landscape C3:C<sup>4</sup> ratio over a much larger area (up to 100s of km<sup>2</sup> ) over a shorter timescale than do paleosol carbonates.

Clades other than Poaceae (grasses) also include lineages that have independently evolved C<sup>4</sup> photosynthesis and may be ecologically dominant, notably Cyperaceae (sedges) (Sage et al., 2011), and some marine and lacustrine isotopic records of C<sup>4</sup> biomass likely reflect local abundance of C<sup>4</sup> plants other than grasses (e.g., Schefuß et al., 2011; Ivory and Russell, 2016). However, the herbaceous component of the spatially extensive modern tropical and temperate grasslands and savannas are overwhelmingly dominated by C<sup>4</sup> grasses (Still et al., 2003; Lehmann et al., 2011). Relatively few studies have compared carbon isotopic data and plant microfossil assemblages from the same terrestrial sites, but those that have (e.g., McInerney et al., 2011; Strömberg and McInerney, 2011; Cotton et al., 2012; Chen et al., 2015; Smiley et al., 2016) suggest that C<sup>4</sup> isotopic signatures in terrestrial sediments are dominated by grasses and not by other clades, such as sedges, that also evolved C<sup>4</sup> photosynthesis. Based on this reasoning, we assume that the terrestrial isotopic records analyzed herein primarily reflect the abundance of C<sup>4</sup> grasses in Neogene ecosystems.

## CARBON ISOTOPE DATA COMPILATIONS

The dataset of paleosol carbonate δ <sup>13</sup>C-values was compiled from publications and datasets in DLF's personal collection and from extensive (but potentially not exhaustive) bibliographic searches using Google Scholar and the GeoRef bibliographic database maintained by the American Geosciences Institute accessed via the University of Minnesota Libraries during 2011–2013. Search terms included combinations of paleosol, pedogenic, carbon<sup>∗</sup> , isotop<sup>∗</sup> , soil, Neogene, Quaternary, Miocene, Pliocene, and Pleistocene, where <sup>∗</sup> is a wildcard operator appropriate for each database. To be included, published datasets had to include specific geographic information for measured sections or individual samples and radiometric ages and/or precise chronostratigraphic assignments in the Miocene (25.0–4.9 Ma) and/or Pliocene (4.9–2.6 Ma) epochs. Most papers included tables of individual δ <sup>13</sup>C values, but for papers in which data were only presented in figures, Data Thief (https://datathief.org/) was used to capture digitally individual data points from figures. Tests of Data Thief on figures with accompanying data tables verified the accuracy of data capture. The paleosol compilation includes data from 44 publications published from 1986 to 2012 (**Data Sheet 1**), including 1,280 values from 189 Miocene sections and 1,379 values from 101 Pliocene sections in Africa, Eurasia (including southeast Europe), North America, and South America (**Figures 1A**, **2**).

The dataset of fossil mammal tooth and ratite eggshell δ <sup>13</sup>C values was initially compiled by Ben Passey (now at University of Michigan, Ann Arbor, MI, U.S.A.), and was completed with the addition of data from publications and datasets in DLF's personal collection and from extensive bibliographic searches using Google Scholar and GeoRef. Search terms included combinations of mammal<sup>∗</sup> , fossil, tooth, enamel, apatite, bioapatite, carbon<sup>∗</sup> , isotop<sup>∗</sup> , Neogene, Quaternary, Miocene, Pliocene, and Pleistocene. Inclusion of papers and data acquisition followed the criteria for paleosol carbonates. The biomineral compilation includes data from 57 papers published from 1994 to 2013 (**Data Sheet 2**), including 1,853 values from 173 Miocene faunal sites and 792 values from 52 Pliocene faunal sites in Africa, Eurasia, North America, and South America (**Figures 1B**, **3**). Paleosol carbonate and biomineral data were also compiled for Pleistocene localities (not exhaustively), and these are summarized in **Figures 2**, **3** but not analyzed here.

To compare published δ <sup>13</sup>C values to the vegetation model results, we first had to correct the δ <sup>13</sup>C values for secular change in the δ <sup>13</sup>C of atmospheric CO<sup>2</sup> over the Miocene and Pliocene (Passey et al., 2009; Tipple et al., 2010), which causes the end-member C3- and C4-values to vary through time. To do this, we binned the data in both compilations into 0.5 Myr intervals (approximating the Miocene-Pliocene boundary as 5.0 Ma and the end of the Pliocene as 2.5 Ma) and used the estimates of the δ <sup>13</sup>C of atmospheric CO<sup>2</sup> (based on the δ <sup>13</sup>C of benthic foraminifera) in those 0.5 Myr bins to correct the measured δ <sup>13</sup>C values to an equivalent value assuming a "modern" δ <sup>13</sup>C value of −8.0% (following Passey et al., 2002, 2009). The modern mean values for C<sup>3</sup> and C<sup>4</sup> plants, the apparent enrichment factors for each photosynthetic pathway, and the uncertainties for the mean values and enrichment factors served as inputs to the linear mixing model IsoError 1.04 (Phillips and Gregg, 2001), which were used to calculate the mean percentage C<sup>3</sup> and C<sup>4</sup> biomass for each corrected δ <sup>13</sup>C value, as well as the standard deviation for each mean value and the 95% confidence interval for each estimate. The enrichment in <sup>13</sup>C by C<sup>3</sup> photosynthesis during fixation of CO<sup>2</sup> decreases with increasing irradiance and water deficit (Cernusak et al., 2013), reducing the spacing between C<sup>3</sup> and C<sup>4</sup> end-members, but given the temporal and geographic scope of the analyses here and for the sake of simplicity, we did not consider the potential influence of light availability and water deficits on δ <sup>13</sup>C value of C<sup>3</sup> plants in our interpretations of the carbon isotopic data.

## PALEOCLIMATE SIMULATIONS

We forced the vegetation models with simulated monthly climates from the HadCM3 family of coupled OAGCMs (Valdes et al., 2017) for Miocene (HADCM3BL-M2.1; Bradshaw et al., 2012) and Pliocene (HADCM3B-M1; Bragg et al., 2012; Conn, 2012, unpublished Master's thesis) paleogeography. Details of these models are given in Valdes et al. (2017). Grid cells are 2.5◦ latitude by 3.75◦ longitude and outputs taken from the OAGCM for each grid cell include monthly temperature, precipitation, and humidity. We used two Miocene simulations: a low pCO<sup>2</sup> case (280 ppmV), which is close to the best fit for proxy data during the late Miocene (Foster et al., 2017) and a high pCO<sup>2</sup> case (401 ppmV), which could represent conditions just after the Miocene climatic optimum (ca. 16.5–15 Ma) (e.g., Kürschner et al., 2008). We used three Pliocene simulations: a low pCO<sup>2</sup> case (280 ppmV), a moderate pCO<sup>2</sup> case (405 ppmV, considered here equivalent to the Miocene high pCO<sup>2</sup> case), and a high pCO<sup>2</sup> case (560 ppmV) that represents high pCO<sup>2</sup> forcing of the Mid Pliocene Warm Period (ca. 3.3–3.0 Ma) (e.g., Haywood et al., 2016). The SDGVM averages the climates for neighboring grid cells and restricts humidity to the range 30–95%, and for consistency these averaged climates were used for all three approaches.

## VEGETATION MODELS

We used the simulated climates to predict dominance or abundance of C<sup>4</sup> grasses in each grid cell employing three increasingly complex vegetation models. The simplest model is the crossover temperature model. A hallmark of C<sup>4</sup> plants is their dominance in high-light and high-temperature environments such as grasslands and savannas (Long, 1999; Sage et al., 1999). The dominant process-based explanation for environmental controls on C<sup>3</sup> and C<sup>4</sup> grass distributions is the crossover

temperature hypothesis, which is based on the different quantum yields of grass species using each pathway (Ehleringer, 1978; Ehleringer et al., 1997). The quantum yield describes how much CO<sup>2</sup> is absorbed by a leaf compared to the amount of APAR it receives during light limited photosynthesis (Ehleringer and Björkman, 1977; Ehleringer, 1978; Collatz et al., 1998). The quantum yield in C<sup>3</sup> plants decreases with increasing leaf temperature, and increases with CO<sup>2</sup> at a given leaf temperature, essentially reflecting the influence of these factors on photorespiration (Ehleringer and Björkman, 1977; Pearcy and Ehleringer, 1983; Collatz et al., 1998; Sage, 2004). By comparison, the quantum yield of C<sup>4</sup> plants is relatively constant across a range of temperatures and CO<sup>2</sup> levels due to the C<sup>4</sup> carbon-concentrating mechanism. The point at which the quantum yield of C<sup>3</sup> grasses equals the quantum yield of C<sup>4</sup> grasses was defined as the "crossover temperature" (Ehleringer et al., 1997; Collatz et al., 1998). At temperatures below this cutoff, C<sup>3</sup> grasses should have a higher capacity to fix carbon, while C<sup>4</sup> grasses should have higher capacities at temperatures above this cutoff. The crossover temperature model simplifies physiological differences between C<sup>3</sup> and C<sup>4</sup> plants but has nevertheless been fairly successful in predicting large-scale distributions (Collatz et al., 1998; Still et al., 2003; Griffith et al., 2015).

For each of the five paleoclimate simulations, we estimated crossover temperatures based on the assumed values of atmospheric pCO2. Representative photosynthetic fluxes were predicted using the coupled C<sup>3</sup> and C<sup>4</sup> leaf photosynthesis and stomatal conductance models of Collatz et al. (1991, 1992). Parameter values, such as maximum carboxylation rates (Vmax) and temperature response functions, were taken from Sellers et al. (1996). Vmax for C<sup>3</sup> grasses was assumed to be 90 <sup>µ</sup>mol m−<sup>2</sup> s −1 at 298 K, and 30 µmol m−<sup>2</sup> s −1 for C<sup>4</sup> grasses at 298 K. These models, which are described in detail elsewhere, estimate gross photosynthetic rates as a function of temperature, relative humidity, insolation, and the partial pressures of atmospheric CO<sup>2</sup> and O2. Crossover temperature values were estimated for light-limited conditions and represent the leaf temperature where C<sup>3</sup> and C<sup>4</sup> photosynthetic rates were equal. Crossover temperatures calculated for the 280, 401/405, and 560 ppmV cases are 14, 20, and 24◦C, respectively.

For each grid cell in each paleoclimate simulation, the crossover temperature model as used here is simply the number

FIGURE 2 | Summary of the compilation of paleosol carbonate δ <sup>13</sup>C values. Number of localities (A) and δ <sup>13</sup>C values (B) in each 0.5 Myr bin. All δ <sup>13</sup>C values by mid-point of 0.5 Myr age bin (C) and by latitude of sample site (D). Mean and standard deviation of estimated percent C4 biomass for each locality by midpoint of 0.5 Myr age bin (E) and by latitude of site (F).

of months for which the temperature is equal to or greater than the calculated crossover temperature (based on the assumed pCO<sup>2</sup> in each OAGCM simulation), which is an estimate of the number of months that production of C<sup>4</sup> grasses should be favored, all else being equal. Our use of the crossover temperature concept is intentionally simplistic as a first step in estimating the dominance or abundance of C<sup>4</sup> grasses.

Our second vegetation model, the Collatz model (Collatz et al., 1998), is slightly more complex, and includes a more reasonable modern climate threshold for C<sup>4</sup> grass dominance that accounts for both crossover temperature (22◦C for contemporaneous pCO2) and the moisture necessary for plant growth (≥25 mm of precipitation in a month that meets the crossover temperature criterion). Thus, for each grid cell in each paleoclimate simulation, the Collatz model determines the number of months for which the modeled temperature is equal to or greater than the calculated crossover temperature and for which precipitation is also equal to or >25 mm. The inclusion of a precipitation screen effectively limits the months considered to the growing season.

Finally, we used the Sheffield Dynamic Global Vegetation Model (SDGVM; Woodward and Lomas, 2004) to estimate the proportion of C<sup>4</sup> biomass in each OAGCM grid cell. Inputs for the SDGVM are monthly temperature, precipitation, and humidity (taken from each OAGCM simulation) and soil texture and dynamics. The SDGVM, which has been used extensively to model both modern and past vegetation (e.g., Beerling and Woodward, 2001; Bond et al., 2005; Scheiter et al., 2012), simulates photosynthesis, respiration, and transpiration for six natural PFTs, and soil carbon and nitrogen dynamics. PFTs include broadleaf deciduous trees, needleleaf deciduous trees, broadleaf evergreen trees, needleleaf evergreen trees, C<sup>3</sup> grasses, and C<sup>4</sup> grasses. Annual C<sup>3</sup> and C<sup>4</sup> grasses are the primary successional PFTs. If the climate permits, the four woody PFTs may acquire land cover at the expense of grasses over decadal and longer timescales, but trees are killed by several processes, including fire. At the beginning of each model year, the probability of fire determines the fraction of land which will be available for grass encroachment at any gridpoint. The SDGVM then assigns cover to C<sup>3</sup> and C<sup>4</sup> grasses based on their relative NPP, without explicit reference to a crossover temperature (Nemani and Running, 1996). With regard to this study, limitations of the SDGVM include the lack of savannaspecific PFTs and a detailed treatment of fire and herbivory. Although NPP and leaf area index are calculated daily and outputted monthly in the model, biomass is updated only at the end of the year, restricting our ability to examine seasonal root and stem growth. Here, we use average annual leaf area index for each PFT, which scales with biomass, to estimate C<sup>4</sup> grass biomass as a percentage of total simulated plant biomass in each grid cell. This third vegetation model represents a considerable increase in complexity relative to the crossover temperature and Collatz models as used here.

To compare the δ <sup>13</sup>C data to the model estimates of C<sup>4</sup> dominance or abundance, we used statistical summaries of all Miocene and Pliocene δ <sup>13</sup>C values for fossil biominerals and for paleosol carbonates separately in each grid cell. Given the global distribution of Miocene and Pliocene continental rock units and associated terrestrial vertebrate faunas, most continental grid cells in the model domains did not have carbon isotope records from Miocene or Pliocene sites, so by necessity only a subset of model grid cells are analyzed in each time interval. When aggregated at the geographic scale of grid cells, the mean and median δ <sup>13</sup>C values are virtually indistinguishable and yield similar results, so below we present only results for mean values. Because the observed and modeled results are not necessarily normally distributed, we used the non-parametric Kendall rank correlation coefficient (Kendall's tau, τ ) to test for statistically significant correlations between % C<sup>4</sup> biomass based on mean δ <sup>13</sup>C values in grid cells and the modeled estimates of C<sup>4</sup> dominance (number of months that C<sup>4</sup> is favored) and biomass (SDGVM output). Kendall's τ ranges from +1 for perfect positive correlation of rank orders of the data to −1 for perfect inverse correlation of ranks, and a τ of 0 indicates complete lack of correlation of ranks. Kendall's τ was calculated using the Stats package in R (R Development Core Team, 2013) and in JMP Pro 13.

## RESULTS

Comparison of the percentage of C<sup>4</sup> biomass estimated from paleosol carbonate δ <sup>13</sup>C values with the number of months favoring C<sup>4</sup> photosynthesis based on the crossover temperature model is presented in **Figure 4**. Kendall's rank correlation coefficients (τ ) for the two Miocene models are positive and low, but the correlation for the 401 ppmV case is statistically significant (p = 0.016). Kendall's τ for the Pliocene models are all statistically significant, and the correlations are stronger than for the Miocene models, but they are all negative, indicating that grid cells with more months above the modeled crossover temperature have generally lower δ <sup>13</sup>C values for paleosol carbonates, suggesting less C<sup>4</sup> biomass. Comparison of the results for the biomineral data with the crossover temperature model results are presented in **Figure 5**. Kendall's τ values are positive for all models, and most are statistically significant. The rank correlations are generally stronger than for the paleosol carbonate results, with τ values up to 0.395.

Comparison of the paleosol carbonate data with the Collatz model are presented in **Figure 6**. The results are similar to those for the crossover temperature model. The Miocene cases have positive τ values that are statistically significant, but the Pliocene cases have statistically significant, negative correlations between percentage of C<sup>4</sup> biomass estimated from the paleosol carbonate δ <sup>13</sup>C values and the number of months that favor C<sup>4</sup> photosynthesis. Comparison of the results for the biomineral data with the Collatz model results are presented in **Figure 7** and are generally similar to the comparison of the biomineral data with the crossover temperature model (**Figure 5**): all τ values are positive and several are statistically significant, but the rank correlations are generally weak with most having τ values of 0.062 to 0.200.

Comparisons of the paleosol carbonate data with the SDGVM output for the two Miocene and three Pliocene cases are presented in **Tables 1**, **2**, respectively. For both Miocene cases,

all correlation coefficients are positive and many are statistically significant, although generally the correlations are weak (τ ranges from 0.052 to 0.322 for the 280 ppmV case and from 0.053 to 0.264 for the 401 ppmV case). For both sets of results, the strongest correlations are for the percentage of total annual biomass that is C<sup>4</sup> grasses (highest τ for both models) and percent C<sup>4</sup> cover. For the Pliocene cases, many correlations are statistically significant, but negative, and the correlations are generally weak.

Comparisons of the biomineral data with the SDGVM output for the two Miocene and three Pliocene cases are presented in **Tables 3**, **4**, respectively. For both Miocene cases, all correlations are positive and many are statistically significant, particularly for the 401 ppmV case. The correlations are generally weak, as for the paleosol carbonate data, and the maximum correlation coefficients are slightly lower than for the Miocene paleosol carbonates (τ of ca. 0.22–0.27 for percentage C<sup>4</sup> grass biomass and cover). The rank correlations of the biomineral estimates of C<sup>4</sup> biomass for the three Pliocene cases are the strongest of any of our comparisons, and most are statistically significant. The correlations are negative for three of the SDGVM output parameters for the 280 ppmV case, and otherwise all correlations are positive. The strongest rank correlations are moderately strong (τ ca. 0.4 for all three cases for percentage C<sup>4</sup> grass biomass and cover), particularly compared to most of the other results. It is also noteworthy that these significant, positive τ values are only for C4-specific model outputs, whereas summary values for all PFTs such as total annual biomass do not show significant positive correlations with proxy data.

## DISCUSSION

Our analyses indicate that none of the models have strong explanatory power for the paleo-proxy data at a global scale, with rank correlations reaching a maximum of ca. 0.4. Thus, the degree of process model complexity does not appear to matter substantially for better hindcasting of C<sup>4</sup> grass dominance, at least not when the data from different regions are analyzed together. In other words, it is not clear which of temperature, length of growing season, and ecological interactions (as variously simulated in the three models) mattered most for Neogene C<sup>4</sup> grass distribution and which other factors may have been equally or more important.

Nevertheless, some general conclusions can be drawn. First, many of the comparisons for which ranks are positively correlated are indeed statistically significant, and these results necessarily imply non-random associations between the measures of C<sup>4</sup> dominance and abundance and the predictions based on the three models. In other words, they provide evidence that climate, as the driver of our vegetation models, played a role in the distribution and abundance of C<sup>4</sup> plants (primarily grasses) during the Neogene. However, most of the positive correlations, including those that are statistically significant, are weak to moderate (all τ -values are <0.44 and most for the

Miocene are <0.20). These results indicate that factors other than climate and atmospheric pCO<sup>2</sup> were as or more important controls on the dominance and abundance of C<sup>4</sup> plants during the Neogene, at least to the extent that our various models capture climate and atmospheric CO<sup>2</sup> impacts on vegetation physiology and dynamics. That the controls on the distribution of C<sup>4</sup> plants are multifaceted is well-understood. However, climate and atmospheric pCO<sup>2</sup> are the two parameters for which we have the most reliable methods of quantitative estimation in the geological record. Other factors, such as fire frequency and intensity and herbivore disturbance, cannot yet be estimated quantitatively in meaningful ways for most terrestrial sites. Thus, our results highlight the current limits on our ability to make sense of the isotopic proxy record of C<sup>4</sup> abundance in the past using our understanding of modern processes.

Recent studies in modern ecosystems have revealed a broad, but imperfect fit between modeled vegetation output and (proxy) vegetation data in terms of C<sup>4</sup> biomass and dominance (e.g., Ehleringer et al., 1997; Collatz et al., 1998; Cramer et al., 2001; Sitch et al., 2003; Still et al., 2003; Lehmann et al., 2011; Ardö, 2015). Although climate undoubtedly plays a major role in driving the balance between (C4) grasses and trees in low- to mid-latitude ecosystems (e.g., Sankaran et al., 2005; Hirota et al., 2011; Staver et al., 2011), differential responses to climate in the distribution of savannas across continents point to the role of both disturbance (fire, herbivory) and ecosystem history (e.g., biogeography, soil type) in shaping vegetation dynamics (e.g., Bond, 2008; Lehmann et al., 2011). In light of these at least generally consistent results, the agreement between the carbon isotope data compilations and the vegetation model estimates of C<sup>4</sup> dominance and abundance is poor despite the frequency of positive and statistically significant results. For many of the cases, all three models predict strong dominance by, or high abundance of, C<sup>4</sup> grasses in many grid cells for which the isotopic data suggest only modest C<sup>4</sup> biomass. Additionally, in many cases the models predict essentially no C<sup>4</sup> biomass for many grid cells, for which estimates based on isotopic data range from 0% to almost 80% mean C<sup>4</sup> biomass. We used the average end member δ <sup>13</sup>C values for C<sup>3</sup> and C<sup>4</sup> plants to estimate the percentage of C<sup>4</sup> biomass, which could lead to overestimates of C<sup>4</sup> biomass from sites which experienced high irradiance and/or water deficits. Those conditions lead to a net decrease in carbon isotope fractionation during fixation of CO<sup>2</sup> by C<sup>3</sup> plants (Cernusak et al., 2013), which decreases the isotopic offset between C<sup>3</sup> and C<sup>4</sup> plants. However, even if we adjusted some or all values to account for this possible effect, the decrease in the abundance of C<sup>4</sup> biomass would be no >20% and many grid cells modeled as having no C<sup>4</sup> grasses would still have mean observed δ <sup>13</sup>C values consistent with 50–60% C<sup>4</sup> biomass.

Another important observation from our study is that the results for the biomineral data (predominantly from fossil mammalian herbivore tooth enamel) are stronger than for the paleosol data in terms of both frequency of statistical significance, and the strength and direction of the rank correlations. For example, many of the results for all three vegetation models yield negative correlations with the estimates of C<sup>4</sup> biomass from



Bold indicates p-values that are statistically significant at α = 0.05.

the paleosol carbon isotope data, and many of these cases are statistically significant, which is the opposite of expectations. Most striking is the uniformity of negative correlations for all of the comparisons with the Pliocene paleosol carbonate data, including for the most sophisticated of our vegetation models, the SDGVM. One explanation for the stronger correlations with biomineral isotopic data is the different biases inherent in the two types of proxy data we analyzed. The isotopic signal in vertebrate herbivores integrates vegetation growing across the landscape, whereas paleosols are essentially time-integrated point samples (Griffith et al., 2017). It may also be that consumption of aboveground biomass by herbivores better mirrors the variation in total C<sup>4</sup> biomass than paleosol carbonates, the carbon for which derives mostly from soil- and root-respired CO2, thus more closely reflecting underground biomass, as in modern soils (Angelo and Pau, 2016). However, the fact that only the Pliocene data show contrasting correlations for the different types of isotopic data suggests that broad formational factors are not sufficient as an explanation.

## Sources of Uncertainty in Model-proxy Data Comparisons

Despite these weak and partly contradictory results, we see no reason to question the utility of carbon isotopic data for inferring the abundance of C<sup>4</sup> plants in soil biomass and aspects of the diet of primary consumers (see discussion in Cerling and Quade, 1993; Kohn and Cerling, 2002; Cerling et al., 2010). However, interpreting either type of isotopic data in terms of landscape scale patterns of C<sup>4</sup> abundance is perhaps not as straightforward as is often assumed. As discussed above, isotopic data from paleosol carbonates are likely to be much less spatially averaged (and more temporally averaged) than isotopic data from vertebrate teeth. Similarly, isotopic data from herbivore teeth suffer from a number of unique biases. For example, carbon isotope data for consumers for a single locality are subject to many filters that blur or obscure how the distribution of diets in terms of C<sup>4</sup> consumption by individuals and averages for species relate to the abundance of C<sup>4</sup> grasses locally or regionally. First, herbivores may not sample the vegetation evenly or completely, so that even if the fauna could be analyzed in its entirety, it would not provide an unbiased isotopic picture. Second, taphonomic factors can alter the taxonomic composition of faunas, both in terms of presence or absence of species in the fossil record and also relative abundance and evenness of species relative to the original living faunal community (Behrensmeyer et al., 2000). Sampling of fossil taxa rarely accounts for the relative abundance of species in the fossil assemblage, which could be an important consideration for interpreting distributions of δ <sup>13</sup>C values or species or faunal mean values in relation to the inferred nature of the ancient habitat. Third, in modern ecosystems in which essentially all grasses are C<sup>4</sup> and all shrubs and trees are C3, such as much of low elevation east African savannas today, using δ <sup>13</sup>C values to classify consumers into distinct dietary categories (C<sup>4</sup> grazer, mixed feeder, C<sup>3</sup> browser, closed canopy C<sup>3</sup> browser) and then to classify faunas into proportions of species in each category, as Cerling et al. (2015) did, is a possible way to convert consumer isotopic data into quantitative models of biome type. However, in regions that have ecologically meaningful abundance of C<sup>3</sup> grasses, this approach may not work as well because some grazers will have a C<sup>3</sup> signal. These filters and biases mean that, for application to fossil faunas, careful consideration of feeding ecology, taphonomy, and sampling are necessary. Nevertheless, the isotopic proxies for C<sup>4</sup> are reliable in principle, if somewhat more nuanced than they are commonly treated, so we do not think that the isotopic data themselves are necessarily problematic.

Similarly, we do not question the basic concepts of the crossover temperature or the modern climatic threshold proposed by Collatz et al. (1998) as predictors of the physiological advantage of C<sup>4</sup> plants under the right environmental conditions. The fact that neither our simple crossover temperature model nor the only slightly more complex Collatz model performs well suggests that pCO2, temperature, and moisture availability alone are not adequate to explain the history of C<sup>4</sup> grasses, despite the fact that the Collatz model does reasonably well at predicting the modern distribution of C3, C4, and mixed grasslands today (Collatz et al., 1998). Other local to regional biotic and abiotic factors on each continent must also have been important—and continue to be important—as has been suggested before (Fox and Koch, 2003; Edwards et al., 2010; Lehmann et al., 2014). Even the SDGVM results are not particularly compelling, despite TABLE 2 | Kendall's rank correlation coefficients (τ ) and p-values for comparisons of estimates of mean percent C4 biomass from δ <sup>13</sup>C-values of Pliocene paleosol carbonate with results of three cases of the SDGVM.


Bold indicates p-values that are statistically significant at α = 0.05.

the attempt in this model to account for a broader range of environmental variables and a more nuanced set of plant life strategies, and a demonstrated ability to reasonably recreate modern vegetation patterns (Cramer et al., 2001; Woodward et al., 2001).

Instead, we argue that an important reason for the relatively weak correlations between model output and proxy-data could be the difference in temporal and spatial scale and resolution between the paleosol and consumer records on the one hand and the OAGCM grid cells on the other. Individual paleosol samples are essentially point samples and at the scale of the OAGCM grid cells even multiple stratigraphic series of samples from a field area are still essentially point samples. The few cases of extensive lateral sampling at the outcrop and field area scale have shown considerable heterogeneity in the abundance of C<sup>4</sup> grasses on length scales of 10s to 1,000s of meters that reflect the influence of landscape processes, even in the context of regional transitions from C<sup>3</sup> to C<sup>4</sup> dominated ecosystems (e.g., Behrensmeyer et al., 2007 for the classic Siwaliks record of Pakistan). Landscape position, local topography, and edaphic factors such as grain size and clay content are critical controls on soil drainage and the potential for carbonate precipitation in addition to climate (Birkeland, 1999), so the paleosol record may not be representative of the average conditions across entire OAGCM grid cells. Vertebrate consumers collectively sample the landscape much more broadly than do soils, but only some of the largest species might have had individual home ranges that approach the scale of the model grid cells and we presume that relatively few species had long distance migration routes, as is the case for most modern mammals (Harris et al., 2009). In addition, as mentioned above, even these large herbivores are unlikely to have sampled the landscape in a perfectly representative way. Given that the consumer data reflect a mix of the feeding preferences of individuals and the typical behavior of populations and species, it is therefore possible that the biomineral dataset also does not completely reflect average conditions at the scale of every OAGCM cell (although our strongest results are for the Pliocene biomineral data in relation to the SDGVM output). Thus, the proxies may be sampling small parts of the area of a grid cell relatively accurately, but those may not be representative of the entire area of the grid cell, while the SDGVM results may also be relatively accurately predicting average conditions, but at the scale of the grid cell. Spatial averaging of nonlinear physiological responses to temperature and CO<sup>2</sup> ("Jensen's inequality"; Denny, 2017), would result in over- or underestimation of C<sup>4</sup> dominance and abundance in the large grid cells of the SDGVM.

The temporal scales of our comparisons may also not be commensurate. The isotopic datasets both span millions of years and are unevenly distributed in time and space. Moreover, in some areas (e.g., south Asia), the isotopic record in either the Miocene or the Pliocene or both has known patterns of temporal


Bold indicates p-values that are statistically significant at α = 0.05.

variability. The OAGCM outputs for the late Miocene and the Pliocene are single time slices that necessarily represent climate over millions of years for each interval and use geographic boundary conditions in the Miocene and Pliocene with their own uncertainties. Furthermore, global and, more importantly, regional climates were likely more dynamic during both periods in ways that could have been important for the evolution of C<sup>4</sup> grasslands. For the Pliocene, the multiple models at different pCO<sup>2</sup> may capture some aspects of the global scale dynamic of the transition into and out of the Mid Pliocene Warm Period (see Haywood et al., 2016). Also, given the relative brevity of the Pliocene compared to the Miocene, the geographic boundary conditions during the Pliocene may not have changed enough to matter. However, that is probably not the case for the Miocene given the histories of active tectonism in western North America, east Africa, and south Asia during this time, and the known influence of changes in orography on climate. These considerations highlight the fact that the paleoclimate simulations that underlie our vegetation models may be inaccurate in specific grid cells or even over regions in ways that contribute to the relatively low correlations between data and models that we observe.

Similarly, despite the comparative sophistication of the SDGVM, various factors may complicate its ability to model ancient ecosystems in a ways that limit the comparisons we have made here, particularly with regard to how grass-dominated ecosystems are modeled. The SDGVM does not have savannaspecific FTs such as fire-adapted, shade-intolerant trees, or a distinct shrub FT. The treatment of fire in this implementation of SDGVM is likely another important factor, given the prominence of fire as a control on modern grass-dominated ecosystems (e.g., Bond et al., 2005). Fire probability is calculated at the beginning of the year using climate envelopes and fuel loads depend on biomass, which is updated annually, so the model is relatively insensitive to vegetation dynamics; in addition, fire intensity is not modeled. The PFTs do not differ in flammability and have no fire-specific adaptations such as thick bark. The treatment of fire may be a reason in some cases for the SDGVM underpredicting C<sup>4</sup> biomass in areas that the isotopic data indicate abundant C<sup>4</sup> biomass. The model does not include herbivory, which like fire is a vital factor in maintaining grassy biomes, and has been implicated in the origins of savannas in Africa (Charles-Dominique et al., 2016). In the implementation we used, climate envelopes directly control tree distribution but only indirectly control grasses via their net primary production. The SDGVM does not have competition between the PFTs, except for crude competition for water resources, but without hydraulic redistribution. These points are not intended as criticisms of the SDGVM, since all models make simplifications. However, these factors may limit how well the SDGVM can model the types of biomes that were dominated by C<sup>4</sup> grasses today and in the past. This is particularly true if a non-analog combination of factors (biotic and abiotic) promoted the spread of C4-dominated vegetation during the Neogene.

Finally, we note that all three vegetation models are underlain by the simulated climates, which have uncertainties that may be equally high or even higher than those of the SDGVM. Indeed, the climate simulated can be highly dependent on the geographical boundary conditions assumed (e.g., Zhang et al., 2012; Brierley and Fedorov, 2016; Ahlström et al., 2017). Validation of the climates simulated by the OAGCM are complicated by the same limitations of the paleoclimate proxy record that preclude quantitative analysis of the abiotic factors that controlled the history of C<sup>4</sup> grasslands globally over the Neogene. However, numerous model-proxy data comparisons for the Cenozoic attest to the inability of OAGCMs to reproduce climatic events in the past (e.g., Herold et al., 2010; Lunt et al., 2012), with the Paleocene-Eocene Thermal Maximum as a particularly well-studied example (McInerney and Wing, 2011). These model-data mismatches indicate fundamental gaps in our understanding of Earth's climate system that have to be filled before more meaningful comparisons can be made between simulations of past climates and proxy data for past vegetation.

## CONCLUSIONS

In the end, given all of the possible complications in our data-model comparisons, as well as the complexity of the regional and local processes involved in the evolution of C<sup>4</sup> grasslands over the Neogene (and Quaternary), perhaps the statistically significant comparisons with rank correlation TABLE 4 | Kendall's rank correlation coefficients (τ ) and p-values for comparisons of estimates of mean percent C4 biomass from δ <sup>13</sup>C-values of Pliocene biominerals with results of three cases of the SDGVM.


Bold indicates p-values that are statistically significant at α = 0.05.

coefficients in the 0.2–0.4 range should actually be viewed as good agreement with the models, particularly the SDGVM. The models are intentional simplifications and abstractions of complex interaction properties and processes, and the historical patterns in the carbon isotopic (and other paleontological and geological) records are complicated in both space and time. Our results are consistent with prior findings that global changes in pCO<sup>2</sup> and global to regional changes in temperature and precipitation alone are not sufficient to explain the histories of C4-dominated ecosystems since the beginning of the Neogene. Our results also highlight the need for means to estimate in the geological record other factors that are known from modern studies to be important controls on C<sup>4</sup> abundance, such as fire intensity and herbivore disturbance. Fire and herbivory should lead to model under-predictions of C<sup>4</sup> abundance, whereas relying only on climate and CO<sup>2</sup> might lead to over-predictions of C<sup>4</sup> abundance. Since abundance of C<sup>4</sup> grasses depends on disturbance factors that limit tree growth, then model omission or over simplistic representation of herbivory and fire can lead to an underestimation of C<sup>4</sup> grass dominance. Similarly, if C<sup>4</sup> grass dominance is simulated solely on the basis of physiological differences with C<sup>3</sup> plants (the crossover and Collatz models) without considering the ecological interactions between grasses and trees, then models will tend to overestimate C<sup>4</sup> grass abundance. Those aspects of the SDGVM that may complicate

its ability to predict grassy biomes suggest that additional model complexity may also be necessary, and regional scale models with higher spatial resolution for multiple time steps through the Neogene and Quaternary are needed to address the scale issues we identified here. However, increasing model complexity and resolution may not yield better model-data comparisons if local ecological interactions and evolutionary history have equally strong effects as abiotic conditions, so that the physiological responses to environment are only one factor among several. Either way, our results should add caution to the reliance on current models for predicting the future of C<sup>4</sup> grasslands, and point to the importance of considering scale when comparing different sources of both ecological and paleoecological data.

## AUTHOR CONTRIBUTIONS

DF, SP, CS, CO, and CJS designed the project. DF compiled the isotopic database, performed the statistical analyses, and wrote the initial draft of the paper. CS calculated the crossover temperatures and Collatz model values. CB and SC performed the climate simulations. LT and DB performed the SDGVM simulations. DF, CS, SP, CO, LT, and CS edited and revised the paper.

## FUNDING

A grant from the National Evolutionary Synthesis Center (NESCent; funded by NSF EF0905606) to CO, CAES, and CJS funded the Origins of C<sup>4</sup> grasslands: a new synthesis of phylogeny, ecology and paleobiology working group that led to this work.

## ACKNOWLEDGMENTS

We thank Ben Passey for sharing his initial work on the biomineral isotope data compilation. We also thank the other members of the NESCent C<sup>4</sup> Grasslands Working Group for

## REFERENCES


valuable discussions that helped shaped this project: T. Michael Anderson, William Bond, Erika Edwards, Elisabeth Forrestel, William Hoffmann, Ben Passey, Nicolas Salamin, and Melinda Smith.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2018.00147/full#supplementary-material

Data Sheet 1 | Palesol carbonate stable isotope data compilation.

Data Sheet 2 | Biomineral stable isotope data compilation.

Basin, Kenya from 4 to 1 Ma. Proc. Natl. Acad. Sci. U.S.A. 112, 11467–11472. doi: 10.1073/pnas.1513075112


Global Productivity, eds J. Roy, B. Saugier, and H. A. Mooney (San Diego, CA: Academic Press), 521–541.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Fox, Pau, Taylor, Strömberg, Osborne, Bradshaw, Conn, Beerling and Still. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Grass Functional Traits Differentiate Forest and Savanna in the Madagascar Central Highlands

Cédrique L. Solofondranohatra1,2 \*, Maria S. Vorontsova<sup>3</sup> , Jan Hackel 3,4 , Guillaume Besnard<sup>4</sup> , Stuart Cable2,5, Jenny Williams <sup>6</sup> , Vololoniaina Jeannoda<sup>1</sup> and Caroline E. R. Lehmann7,8,9

<sup>1</sup> Laboratoire de Botanique, Département de Biologie et Ecologie Végétales, Faculté des Sciences, Université d'Antananarivo, Antananarivo, Madagascar, <sup>2</sup> Kew Madagascar Conservation Centre, Antananarivo, Madagascar, <sup>3</sup> Comparative Plant and Fungal Biology, Royal Botanic Gardens, London, United Kingdom, <sup>4</sup> Laboratoire Evolution and Diversité Biologique, CNRS/ENSFEA/IRD/Université Toulouse III, Université Paul Sabatier, Toulouse, France, <sup>5</sup> Conservation Science, Royal Botanic Gardens, London, United Kingdom, <sup>6</sup> Biodiversity Informatics and Spatial Analysis, Royal Botanic Gardens, London, United Kingdom, <sup>7</sup> School of GeoSciences, The University of Edinburgh, Edinburgh, United Kingdom, <sup>8</sup> Centre for African Ecology, School of Animal and Plant Sciences, University of Witwatersrand, Johannesburg, South Africa, <sup>9</sup> Royal Botanic Garden Edinburgh, Edinburgh, United Kingdom

### Edited by:

Daniel M. Griffith, Oregon State University, United States

#### Reviewed by:

Zehao Shen, Peking University, China Deusdedith M. Rugemalila, Wake Forest University, United States

> \*Correspondence: Cédrique L. Solofondranohatra lovacedrique@gmail.com

#### Specialty section:

This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution

Received: 11 January 2018 Accepted: 23 October 2018 Published: 15 November 2018

#### Citation:

Solofondranohatra CL, Vorontsova MS, Hackel J, Besnard G, Cable S, Williams J, Jeannoda V and Lehmann CER (2018) Grass Functional Traits Differentiate Forest and Savanna in the Madagascar Central Highlands. Front. Ecol. Evol. 6:184. doi: 10.3389/fevo.2018.00184 Grassland, woodland, and forest are three key vegetation types that co-occur across the central highlands of Madagascar, where the woodland has historically been considered as degraded forest. Here, we use grass functional traits to inform our understanding of the biogeography of Malagasy vegetation and the differentiation of vegetation types in the region. We sampled grass community composition at 56 sites across the central highlands of Madagascar spanning grassland, woodland, and forest. We selected seven functional traits known to correlate with different aspects of life history collated via GrassBase (habit, culm type, physiology, leaf consistency, plant height, leaf width, and leaf length) for the 71 constituent species. Via analyses of the beta diversity, rank abundance, functional dispersion, functional group richness, and community phylogenetic structure of grassland communities, we differentiate these vegetation types using plant functional traits. Grassland and woodland are highly similar in grass species composition and dominated by the same species (Loudetia simplex, Trachypogon spicatus, and Schizachyrium sanguineum). In contrast, forest grass species composition significantly differs from both grassland and woodland. Consistent with these species composition patterns, the vegetation types can be distinguished based on physiology, culm type, and leaf consistency. Forests harbor primarily C<sup>3</sup> grasses, which are almost invariably laterally spreading with herbaceous leaves. In contrast, both grassland and woodland species are predominantly tall, caespitose C<sup>4</sup> grasses with coriaceous leaves. Forest grasses are phylogenetically clustered and less diverse than the grassland and woodland communities. Further, we sampled bark thickness of the common woody species occurring in the woodland and forest of the region and found that the relative bark thickness of the woodland tree species was significantly greater than that of forest species from the same genus. We found that the functional traits and architecture of grasses diverge strongly between forest and the grassland-woodland mosaic. We conclude that the woodlands, primarily dominated by Uapaca bojeri Baill., are a form of savanna and not forest as has been previously suggested.

Keywords: functional traits, species diversity, grassland, Poaceae, woodland, forest, savanna, Uapaca bojeri

## INTRODUCTION

Forest and savanna are vegetation types that can co-occur within landscapes and found in areas with similar climate and soils (Bond, 2005; Bond and Parr, 2010; Charles-Dominique et al., 2015). While savanna and forest can appear superficially similar in terms of the density and number of trees, the functional traits of the trees of each vegetation type differ (Ratnam et al., 2011; Griffith et al., 2017). Savanna and tropical grasslands are open canopy environments with variable tree cover and a continuous grassy ground layer (Bond et al., 2003a; Bond, 2008; Lehmann et al., 2009; Bond and Parr, 2010; Edwards et al., 2010). Defining the functional nature of a vegetation type is often linked to a broader issue of vegetation classifications based on structural associations and even climatic definitions (Schimper and Fisher, 1903; Whittaker, 1975; Scholes and Hall, 1996; Woodward et al., 2004), and often almost entirely based on the woody components of vegetation. The Food and Agriculture Organization of the United Nations (FAO) defines forest as land spanning more than 0.5 ha with more than 10% tree canopy cover and trees higher than 5 m (or having the potential to reach a height of 5 m; FAO, 2006). Alternatively, the United Framework Convention on Climate Change (UNFCCC) defines thresholds of 0.5–1.0 ha forested land, 10–30% tree canopy cover, and 2–5 m tree height. While savannas are defined as mixed tree–grass systems with a discontinuous tree canopy in a continuous grass layer (Scholes and Archer, 1997; House et al., 2003), the notion of the limits of tropical grassy vegetation types as being defined by dominance of C<sup>4</sup> grasses in the ground layer is more recent (Bond, 2008). House and Hall (2001) described savannas as a continuum between tropical forests and grasslands and use the term to include all plant communities from treeless grasslands to closed canopy woodlands. In other literature, grasslands encompass non-woody and woody grasslands (savannas, woodlands, shrublands, and tundra) together (White et al., 2000) because they are defined as terrestrial ecosystems dominated by herbaceous and shrub vegetation, maintained by fire, grazing, drought, and/or freezing temperatures.

Persistent confusion around the delineation of vegetation types often centers on the limited ecological information in structural definitions (Ratnam et al., 2011; Parr et al., 2014; Lehmann and Parr, 2016). While structural definitions offer precision, they are inadequate when the functional responses of forest and savanna plant species to processes such as fire, herbivory, and drought differ (Lehmann et al., 2011). Several studies have shown functional trait differences between trees of forests and savannas in South America (Hoffmann et al., 2003, 2005, 2012), Australia (Lawes et al., 2013) and Africa (Charles-Dominique et al., 2015), all demonstrating that savanna trees have morphological and physiological traits which make them adapted to the frequent fires that characterize savannas in high rainfall regions (Ratnam et al., 2011; Charles-Dominique et al., 2015). Hoffmann et al. (2003, 2012) showed that bark thickness is a key trait to attain fire-resistance threshold and savanna species have thicker bark than forest species, thereby reducing the risk of stem death during fire.

However, the composition and traits of the grassy ground layer are often neglected yet responsible for the physiognomic structure of vegetation (Bond, 2008), promoting either fire, or herbivory due to grass productivity and architecture (Knapp et al., 1998; Sage and Monson, 1998; Archibald and Hempson, 2016). On the other hand, C<sup>4</sup> grasses are considered absent from forest understorey, which can be dominated by herbaceous life forms (Ratnam et al., 2011). However, some C<sup>4</sup> grasses are known from forest understories [e.g., Setaria palmifolia (J. Koenig) Stapf, Microstegium nudum (Trin.) A. Camus, and Digitaria radicosa (J. Presl) Miq. (Clayton, 1970)]. There are numerous areas of the world where C<sup>3</sup> grasses either dominate or co-dominate with C<sup>4</sup> grasses and these are still considered savanna [e.g., Neotropical savannas (Cabido et al., 1997; Bremond et al., 2012), and tropical savannas in Africa and Australia (Lattanzi, 2010)]. Overall, the focus of studies examining the functional traits of vegetation in distinguishing forest and savanna has been woody rather than grass species (Hoffmann et al., 2003, 2005, 2012; Ratnam et al., 2011; Lawes et al., 2013; Charles-Dominique et al., 2015). While photosynthetic pathway is no doubt of high relevance, this focus also belies important functional variation among grasses of the same photosynthetic pathway (Forrestel et al., 2015; Simpson et al., 2016).

Across Madagascar, a mosaic of grassland covers 65 % of the island and in the central region, it is interspersed with subhumid forest, along with extensive areas of woodland (Moat and Smith, 2007). In this region, woodlands are dominated almost solely by Uapaca bojeri Baill. (commonly called Tapia) in association with a few other woody species that belong to families endemic to Madagascar, especially Sarcolaenaceae (Leptolaena pauciflora Baker, Sarcolaena oblongifolia F. Gérard, Schizolaena microphylla H. Perrier) and Asteropeiaceae (Asteropeia densiflora Baker, A. labatii G.E. Schatz, Lowry & A.-E. Wolf). Malagasy grasslands and savannas have long been considered as anthropogenically derived, the result of forest degradation via clearing, and altered fire and grazing regimes over the last 2,000 years (Perrier de la Bâthie, 1921; Humbert, 1927; Koechlin et al., 1974; Lowry et al., 1997). As such, woodlands of the central highlands have been historically and persistently classified as degraded forests (Koechlin et al., 1974; White, 1983; Faramalala, 1995) and managed as such through suppression of fire (Kull, 2002). Recent studies however support early origins of these extensive grasslands with their numerous endemic grass, forb, and animal species that evolved locally millions of years before human colonization (Bond et al., 2008; Vorontsova et al., 2016; Hackel et al., 2018), which occurred only a few 1,000 years before present (Pierron et al., 2017). Grassland, woodland and forest co-occur in the central highlands of Madagascar with changes in vegetation types occurring over small scales, often just meters where there are no discernible differences in climate and soils. These mosaics offer an opportunity to study the functional and compositional characteristics of the ground layer among these vegetation types to identify the disturbance regimes each is adapted to. Within Poaceae are a broad suite of architectural and morphological functional traits commonly assessed in woody plants, but little explored in grasses. In this study, we examine the composition and functional traits of grass communities among these three vegetation types to quantify differences in community composition and grass functional traits. This allows us to test whether any differences map to vegetation types, in a region where woodlands have long been considered degraded forests but where limits between woodland and grassland appear idiosyncratic. We then assess the phylogenetic structure of each vegetation type as an indicator of its evolutionary history. Finally, as an independent test of vegetation type differences, i.e., savanna vs. forest, we compare the relative bark thickness of woody species, including closely related species of Uapaca that characterize woodland and forest, to determine whether functional traits known to differ between forest and savanna can corroborate findings of the grass communities. We hypothesize that grassland and woodland harbor a high phylogenetic diversity of grass species which have evolved for a long period before human arrival. We predict a functional similarity between grassland and woodland grass community, where light availability and fire regimes are similar. In contrast, forest vegetation, characterized by a low-light environment in the ground layer, will select for a different grasses community with another set of functional traits.

## MATERIALS AND METHODS

## Study Sites

The grass flora was sampled at 56 sites distributed across the central ecoregion of Madagascar, which measures 192,141 km<sup>2</sup> (Humbert, 1955) (**Figure 1** and **Supplementary Table 1**). From North to South, Ibity, Ambatofinandrahana, Itremo, Ambalavao, Andringitra, Isalo, and Ihorombe were the study locations. Rajeriarison and Faramalala (1999) categorized the climate of these regions as subhumid with a mean annual rainfall ranging from 1,000 to 1,500 mm and 5–7 months of dry season. The elevation ranges from 800 to 1,800 meters and soils are primarily ferrallitic, on sandstone and basement gneiss (Moat and Smith, 2007). Fire regime in the sites is every one to 3 years (Alvarado et al., 2018). In our study, we are going to use "grassland" for treeless vegetation, consisting of continuous C<sup>4</sup> grass layer, and "woodland" for woody grasslands. Forest is defined as closedcanopy vegetation, with no C<sup>4</sup> grass layer. We selected sites considered as undisturbed by the absence of invasive species and overgrazing. Each site was in an area of uniform vegetation that was a minimum of 60 m wide. We sampled 24 sites of grassland, 20 sites of tapia (Uapaca bojeri) dominated woodland with herbaceous layer and 12 sites of gallery forest.

## Data Collection

## Grass Species Composition

Based on the sampling method described by Vorontsova et al. (2016), we collected data on grass species composition at 56 sites; data from 17 sites from Vorontsova et al. (2016) that met the above criteria were also included. The method focuses on grass species composition and frequency and consists of four transects. First, all grass species within a center circle plot of one meter diameter were documented. Then, from the center point, four 25-m transects, each following a random direction (based on a compass bearing), were laid out. Along each transect, the same circular plot of one-meter diameter was sampled at five-meter intervals, representing grass species composition over 16.5 m<sup>2</sup> for a single site (see **Supplementary Table 1**).

### Grass Functional Traits

Across these 56 sites, 71 grass species were recorded, and for each species, seven functional traits were collated: (1) habit (annual, perennial); (2) culm type (prostrate, decumbent, ascendant, erect); (3) leaf consistency (herbaceous, coriaceous); (4) culm height; (5) leaf width; (6) leaf length, and (7) photosynthetic pathway. These traits were derived from GrassBase, a global database of 11,369 standardized grass species descriptions compiled from taxonomic literature (Clayton et al., 2006). These data were supplemented from literature (Bosser, 1969; Van Oudtshoorn, 2002) and by data gathering from multiple herbarium specimens held at the herbarium of the Parc Botanique et Zoologique de Tsimbazaza and the Royal Botanic Gardens, Kew, UK. Photosynthetic pathway data were obtained from the checklist of Osborne et al. (2014); for species with an undetermined pathway, it was inferred from their position in known C<sup>3</sup> or C<sup>4</sup> lineages (Hackel et al., 2018).

The choice of these seven traits was based on their role related to light acquisition and flammability. The annual or perennial habit of the species was recorded. In GrassBase, perennials are described as possessing one or more of the following: rhizome or rootstock; a compact mass of dead basal sheaths remaining from previous year's growth; dormant buds in the basal leaf axils for the following year's growth. Annual grass species are commonly considered more abundant than perennials in disturbed habitats, even those with frequent fire (Grime, 1977; Pitelka, 1977). Culm type was categorized into four types: prostrate (lying flat), decumbent (lying along the ground, but with culm apex upright), ascendant (growing upward after an initial horizontal segment) and erect (fully upright) (Beentje, 2010). A prostrate culm type is expected to aid light interception in shaded environments by minimizing self-shading within the plant via lateral spreading (Falster and Westoby, 2003). Erect life forms tend to be species with intravaginal buds [buds between leaf sheaths, not breaking through them to form stolons and rhizomes as in extravaginal buds (Brejda et al., 1989)], that can increase protection of buds from fire, but also confine the architecture of a grass by preventing lateral spreading as by definition these grass species cannot root at the node (Linder et al., 2018), and are therefore generally associated with environments characterized by frequent fire and high light. Depending on the toughness of the leaf blade, we recorded leaf consistency which was categorized into two groups: coriaceous [leathery, tough (Beentje, 2010)], or herbaceous [with the texture of herb, soft and pliable (Beentje, 2010)] leaves.

## Data Analysis

## Environmental Occurrence of Grassland, Woodland, and Forest

The central highlands were delimited following Humbert's (1955) ecoregion map classifying it as "Central ecoregion." The vegetation map developed by Moat and Smith (2007) was used to define the limits of our three vegetation types. Vegetation

units classified as "plateau grassland—wooded grassland mosaic and wooded grassland—bushland," "tapia forest" and "humid forest" in the vegetation map correspond, respectively to our studied vegetation types grassland, woodland, and forest. Any other vegetation units in the vegetation map were omitted from our analysis. The distribution of the three vegetation types across climatic variables (Bio\_1: annual mean temperature and Bio\_12: annual precipitation) obtained from Worldclim Global Climate Data version 2 (Fick and Hijmans, 2017), within the central ecoregion was plotted (see **Supplementary Figures 1, 2**). The variables were downloaded at 30 s spatial resolution. From these, we have generated datasets per climatic variable for our studied ecoregion. Data for each environmental variable were ordered from lowest to highest and binned in the relevant 1◦C and 50-mm intervals, respectively. Due to the very high density of points within each bin, we randomly subsampled our data into 1,000 points, using 1,000 iterations, without replacement. For the subsampled data within each interval, the probability of the presence of each of the three vegetation types was calculated as the mean of all points within the temperature and precipitation interval and the standard errors were calculated for each proportion. The probability was then plotted against Bio\_1 and Bio\_12 to compare the distribution of grassland, woodland and forest within the region. We fitted our 56 sampled sites within the temperature and precipitation gradient as rugs on each x-axis of the probability plots. For each temperature and precipitation interval, the area of land covered by each vegetation was plotted. From the raw climate data, the area of each vegetation type was obtained by counting the number of points within each interval as each point corresponds to one kilometer square of land.

### Grass Species and Trait Composition

We examined the composition of grass communities from three distinct vegetation types as defined by their woody component (or its absence) and its traits: grassland, woodland and forest. All statistical analyses were performed in R 3.0.2 (R Core Team, 2013). We used permutational multivariate analysis of variance (PERMANOVA) using the Bray–Curtis dissimilarity (see **Supplementary Table 2**) and accounted for spatial autocorrelation by including latitude and longitude as predictor variables to assess differences in plant community composition between the three vegetation types, with the R package "vegan" (Oksanen et al., 2016). Grass species and tribes were ranked by frequency per vegetation type with the "BiodiversityR" package (Kindt and Coe, 2005).

To measure species functional diversity in each vegetation type, we computed the functional dispersion index (FDis) which estimates the dispersion of species in a trait space, based on principal coordinate analysis (PCoA) and calculated on a Gower dissimilarity matrix (Gower, 1971; Laliberté and Legendre, 2010) with the "FD" package (Laliberté et al., 2014). We compared the three vegetation types in terms of FDis with an analysis of variance (ANOVA) while accounting for spatial autocorrelation. A functional group richness (FGR) was used to group grass species based on functional traits, using the same package. The resulting dendrogram was plotted against the functional traits. The number of clusters was identified visually on the resulting dendrogram.

Community phylogenetic structure was assessed per site to test how closely related the average pair of grass individuals were in each vegetation type. We used a Bayesian time-calibrated phylogenetic tree of the 71 species sampled (**Supplementary Figure 3**) extracted from a larger analysis of Malagasy grass species (Hackel et al., 2018). Hackel et al. (2018) did not include Andropogon trichozygus Baker, Hyparrhenia variabilis Stapf and Ischaemum polystachyum J. Presl, which were replaced with the closest available species [Andropogon huillensis Rendle, Hyparrhenia schimperi (Hochst. ex A. Rich.) and Ischaemum koleostachys (Steud.) Hack., respectively]. The species richness (Whittaker, 1972) within each vegetation type was calculated. Phylogenetic relatedness of species within grassland, woodland, and forest was measured via phylogenetic diversity (PD) index (Faith, 1992). The differences in species richness and phylogenetic diversity between vegetation types were assessed with a generalized linear model (GLM) with a Poisson distribution and log link function. The net relatedness index (NRI) and nearest taxon index (NTI) (Webb, 2000; Webb et al., 2002) for each site was calculated with the "picante" package (Kembel et al., 2010). The observed phylogenetic relatedness was compared to the expected pattern using the functions ses.mpd and ses.mntd with the "richness" null model, not taking into account species abundances in each plot, and using 999 randomizations.

### Tree Bark Thickness

The bark thickness and stem circumference of 39 individuals of Uapaca bojeri in woodlands, and 27 samples from 16 species in forest (listed in **Supplementary Table 3**) were measured at Isalo National Park and on wood samples from humid forest in the East of Madagascar, available at the forestry department of Antananarivo University, respectively. The tree species from Eastern humid forest are mostly dominant species in our study region (KMCC, 2012). A t-test was used to compare bark thickness between Uapaca bojeri dominated woodland and forest species, and as bark thickness varies with tree size (Werner and Murphy, 2001; Lawes et al., 2011), measures of bark thickness were transformed into relative bark thickness (the ratio of individual bark thickness to stem diameter; Lawes et al., 2013) prior to analysis.

## RESULTS

## Distribution of Grassland, Woodland, and Forest Across Temperature and Precipitation

Across the central ecoregion mean annual temperature ranges from 11 to 27◦C and mean annual rainfall ranges from 764 to 2,660 mm per year. Here, grassland has a greater than 50% probability of occurrence where mean temperatures range from 16 to 22◦C, and a 90% of probability of occurrence where temperature ranges from 11 to 12◦C and 22 to 27◦C. Grassland

Paniceae tribe.

has a >50% probability of occurrence where rainfall ranges from 764 to 864 mm and from 1,160 to 1,760 mm, and a 90% or greater probability of occurrence where rainfall ranges from 864 to 1,160 mm. The probability of forest occurrence is >50% where temperature ranges from 13 to 15◦C and mean rainfall ranges from 1,810 to 2,010 mm and a 90% greater probability of occurrence where rainfall is >2,010 mm.

In total, the sampled woodland covers 4,816 km<sup>2</sup> which represents only 0.95% of the central highlands total land area (**Supplementary Figures 1B, 2B**). The probability of occurrence of woodland is <5% across the rainfall and temperature gradient, corresponding to the limited land area occupied by this vegetation type. However, the greatest probability of occurrence is at 23 and 24◦C. Our studied sites are primarily located in areas with a medium annual mean temperature (20–23◦C). Woodland is only found in very low precipitation areas with a highest probability of occurrence of 25% where rainfall ranges from 760 to 800 mm.

## Grass Diversity and Composition

A total of 71 grass species from 40 genera were recorded and identified using the Poaceae checklist of Itremo Protected Area (Nanjarisoa et al., 2017) and specimens at K, P, and TAN herbaria. We found between 2 and 19 species per plot, with a mean of 8 species. Species richness was higher in woodland and grassland than in forest, with a significant difference between the three TABLE 1 | Comparison of physical environments, species composition, and traits of dominant grass species and trees in grassland, woodland, and forests.


habitats (GLM with Poisson distribution, P < 0.001). Although there is some spatial autocorrelation, the type of vegetation has a significant effect on species composition according to PERMANOVA (F = 14.42, P = 0.001). The vegetation type and the spatial location explain 33 and 9% of the variance, respectively. However, when separate pairwise comparisons are made, grassland and woodland do not differ significantly in composition (F = 1.82, P = 0.06). There are highly significant differences between both the composition of grassland and forest (F = 19.43, P = 0.001) and woodland and forest (F = 23.95, P = 0.001).

The same three species that are most common in grassland are also the dominant species in woodland: Loudetia simplex (Nees) C.E. Hubb., Trachypogon spicatus (L. f.) Kuntze and Schizachyrium sanguineum (Retz.) Alston (**Figures 2A,B**). All three are absent from forest, which is most commonly characterized by Oplismenus compositus (L.) P. Beauv., O. hirtellus (L.) P. Beauv., O. flavicomus Mez, and Arundinella nepalensis Trin (**Figure 2C**).

Rank-frequency curves per grass tribe also support differences in species composition of grassland and woodland in comparison with forest (**Figures 2D–F**). The two Poaceae tribes Andropogoneae and Tristachyideae dominate grassland and woodland, whereas forest species mostly belong to the tribe Paniceae.

## Grass Traits

Beyond the similarities in taxonomic composition, grassland and woodland grass species were also characterized by similar trait communities (**Table 1**). This is evidenced by the functional group richness classification (**Figure 3**) which identified two distinct groups based on the seven functional traits. The two supported clusters demonstrate a close match with vegetation type.

The first cluster of 44 species (called "Non Forest") is primarily composed of species present in grassland and/or woodland (**Figure 3**). Within this group, 89% are C4, 86% are perennial, 91% are erect and 52.3% have coriaceous leaves. The median values of height, leaf width and leaf length of the species in this group are 100, 0.5, and 30 cm, respectively. The key trait that differentiated this "Non-Forest" from the "Forest" cluster was the architecture of culms (**Figure 3**). Grassland and woodland grass culms are primarily erect, leading to a tufted growth form. Their leaves are narrower but also longer and more coriaceous compared to those in the "forest" group which are wide, herbaceous and flaccid.

Of 19 grass species recorded in forest, 84% cluster in the "Forest" group by traits (**Figure 3**). The three forest grass species [Arundinella nepalensis, Lasiorrhachis hildebrandtii (Hack.) Stapf, and Panicum umbellatum (Trin.)] which do not cluster in this group are C<sup>4</sup> and perennial. Species that also occur in grassland and woodland were assigned to the "forest" cluster (7 of 52 species): Panicum perrieri A. Camus, P. subhystrix A. Camus, P. ibitense A. Camus, P. spergulifolium A. Camus, Festuca camusiana St.-Yves, Styppeiochloa hitchcockii (A. Camus) Cope, and Trichanthecium brazzavillense (Franch.) Zuloaga & Morrone. These species are characterized by traits common amongst forest grasses as all species are C3, and all but one have herbaceous leaves that are shorter and wider. The Panicum species are often found in the sub-canopy of the grassland and woodland layer.

Species dispersion in functional trait space differed significantly between the three vegetation types according to ANOVA (P = 0.014). Grassland communities are significantly less diverse functionally compared to forest communities (P = 0.009) and no significant difference was found between grassland and woodland (P = 0.06). Functional diversity is similar between forest and woodland (P = 0.18).

## Diversity and Phylogenetic Structure

Species richness of grass species is significantly different between the three vegetation types (GLM, P = 0.001). Species richness is the highest in woodland with a significant difference compared to grassland (GLM, P = 0.005) and forest (GLM, P < 0.001; **Figure 4A**). Phylogenetic diversity means are equal to 233.1, 272.3, and 140.1 for grassland, woodland, and forest, respectively. Woodland species are significantly more phylogenetically diverse compared to grassland (GLM, P < 0.001) and forest species (GLM, P < 0.001; **Figure 4B**). When compared to a null model, forest grass communities exhibit a pattern of phylogenetic clustering but there is a significant difference with grassland, only for NTI (ANOVA, P = 0.02). No significant difference is found between grassland and woodland (ANOVA, P = 0.145), or woodland and forest (ANOVA, P = 0.205) NTI. Woodland communities are only slightly clustered (NRI = 0.28; NTI = 0.20) and grassland communities are slightly over-dispersed phylogenetically (NRI = −0.01; NTI = −0.08; **Figures 4C,D** and **Supplementary Figure 3**).

## Woodland and Forest Trees Traits

The three vegetation types also differ in terms of woody species composition and traits (**Table 1**). While there are no woody species in grassland, the woodland is dominated by Uapaca bojeri, associated with other woody species, especially Sarcolaena oblongifolia and Asteropeia densiflora. Tapia trees are shorter (8–12 m; Kull, 2002) than humid forest trees. Humid forest canopy height can reach 30–35 m at low to mid-elevations (Moat and Smith, 2007). Moreover, forest trees have wide and dense canopies, creating a low-light understorey. Woodland on the other hand is characterized by a high-light understorey as trees have open crowns which allow light to penetrate.

Uapaca bojeri has significantly thicker bark than the forest species (t = −16.5, P < 0.001; **Figure 5**). Among the forest species, five samples were from the same genus Uapaca (two U. ferruginea Leandri and three U. densifolia). The bark of U. bojeri (13.5%) is at least twice thicker that in forest Uapaca species (6.15% in U. ferruginea; see **Supplementary Table 3**).

## DISCUSSION

The environmental envelopes of the three vegetation types show that grassland and woodland are most likely to co-occur, while forest has a distinct climatic envelope. Forest occurs at temperature ranging from 13 to 15◦C and precipitation ranging from 1,810 to 2,660 mm, whereas grassland and woodland occur at temperature ranging from 18 to 25◦C and precipitation ranging from 765 to 1,610 mm (**Supplementary Figures 1, 2**). An independent assessment based on herbarium records along environmental variables yielded an extremely similar environmental distribution of the Uapaca bojeri dominated woodland (Werkmeister, 2016). Werkmeister (2016) showed that annual temperature and mean annual precipitation are important factors in shaping the woodland's modern distribution, and the overall ranges of the woodland's environmental envelopes (temperature ranging from 15,1 to 25.1◦C and precipitation ranging from 731 to 1,576 mm) coincide with our results of its occurrence probability. Our results on grassland and woodland distribution in Madagascar are similar to the distribution of savannas on the African continent where mean annual rainfall and rainfall seasonality are the key environmental correlates (Lehmann et al., 2011). The fire regime of the central highland grasslands has a return time of one to 3 years (Alvarado et al., 2018) which is in line with observations made for the African continent where the climate is similar (Lehmann et al., 2014). The locations sampled in this study span a precipitation gradient from 781 to 1,412 mm, similar to the rainfall of savannas in Africa, and corresponds to geographic regions where the three vegetation types can locally co-occur as a mosaic within the landscape. Further, in this study, we sampled 13% of the grass flora of Madagascar [71 out of 550 species (Vorontsova et al., 2016)] and these are taxonomically representative of the central highlands. Hence, the divergence we found between the diversity and functional traits of grass species in woodland and grassland in comparison to forest (**Figures 2**, **3**) are representative of community composition and biogeographic patterns of these vegetation types across the region. Our analyses identified two distinct functional groups, that of grassland and woodland species vs. forest species defined by differentiation in architecture, physiology, and leaf traits, all of which are related to light interception, fire tolerance and flammability (**Figure 3**). The ground layer of Uapaca bojeri dominated woodland is characterized by grasses of the tribes Andropogoneae and Tristachyideae (**Supplementary Figure 3**), the same that characterize central Madagascar grasslands. Studies have shown that the tribe Andropogoneae exhibits a suite of functional traits that are advantageous in regularly burned grasslands (Everson et al., 1988; Morgan and Lunt, 1999; Bond et al., 2003b; Spasojevic et al., 2010; Forrestel et al., 2014; Simpson et al., 2016). The majority of grass species found in woodlands

Frontiers in Ecology and Evolution | www.frontiersin.org November 2018 | Volume 6 | Article 184

species with similar traits from forest and non-forest habitat.

and grasslands are tall perennial C<sup>4</sup> grasses, with coriaceous leaves, that burn frequently (**Figure 3**).

Grass species in forest and woodland showed a pattern of higher dispersion in the functional trait space than in grassland. This indicates a greater environmental heterogeneity in these vegetation types, where species exhibit a wide range of traits to adapt to their habitat (Anderson et al., 2006). Several studies have shown that high light intensity and high temperature environments favor C<sup>4</sup> over C<sup>3</sup> photosynthesis (Ehleringer et al., 1997; Osborne and Sack, 2012), and C<sup>4</sup> grasses dominant in tropical savannas and grasslands (Sage and Monson, 1998; Bond, 2008). C<sup>4</sup> grasses which dominate the woodland generally have coriaceous leaves, linked to traits that retard decomposition (Street et al., 1980) and increase flammability such as high leaf C:N ratios (D'Antonio and Vitousek, 1992; Mouillot and Field, 2005; Bond, 2008; Cardoso et al., 2008). These species are also tall grasses, promoting fire and making the vegetation highly flammable (Fensham et al., 2003; Vigilante et al., 2004; Archibald et al., 2005; Simpson et al., 2016) and are also strong light competitors (Olff and Ritchie, 1998). In contrast, grass species in the forest understory are essentially C<sup>3</sup> members of the tribe Paniceae (**Supplementary Figure 3**). Forest understories are characterized by high moisture and low light conditions (Hoffmann et al., 2003, 2012), and forest grasses are largely annual, creeping, with wide leaves that enable them to optimize light interception (Milla and Reich, 2007) and minimize self-shading (Falster and Westoby, 2003). Leaf size has been demonstrated to be related to light acquisition. A significant positive relationship between leaf size and specific leaf area has been reported (Ackerly and Reich, 1999); and high specific leaf area has been shown to allow species to optimize light interception in shady conditions (Milla and Reich, 2007). A study from a woodland community in Australia has shown that small-leaved woody species were associated with more light in the understorey (Bragg and Westoby, 2002). In contrast, high-light environments of grassland and woodland are dominated by perennial tufted grasses with erect culms. With an erect architecture, these species can protect their meristematic tissue from fire damage with intravaginal buds protected within basal leaf sheaths or underground, and tillers

tightly clustered (Linder et al., 2018). Thus, after fire, such species rapidly resprout (Ripley et al., 2015), and proliferate under the disturbance regime that characterizes grassland and woodlands, where humans have limited influence on the fire regimes (Alvarado et al., 2018).

Tapia woodland grass species have higher species richness and phylogenetic diversity compared to grassland and forest species (**Figures 4A,B**). This could be explained by the environmental gradient within the woodland promoted by heterogeneity in tree cover, enabling co-existence of a higher diversity of grass species (Silvertown, 2004) than in grassland or forest. There is an overall greater range of light environments within the woodland. Like grassland, woodland burns frequently, but the forests rarely burn (Ratnam et al., 2011; Charles-Dominique et al., 2015). Our analyses of community phylogenetic structure found evidence of evenness in grassland, slight clustering in woodland and clustering in forest. In forest, ground layer composition must be adapted to light availability (Charles-Dominique et al., 2015), potentially leading to phylogenetically conserved traits at local spatial scales. Here, most grass species found in forest belong to the "Madagascar shade clade" in the tribe Paniceae (Hackel et al., 2018) and share similar physiological and biological traits (creeping C<sup>3</sup> species with herbaceous leaves). In contrast, it has been suggested that competitive exclusion plays a structuring role if co-occurring species are less closely related than expected by chance (Weiher and Keddy, 1995; Webb, 2000). Grassland and woodland grass species are evenly distributed across the phylogeny, and these communities are characterized by competition for resources and regular disturbance (Bond, 2008). No significant difference was found between the phylogenetic distance to the nearest taxon for each taxon in the sample (NTI) between grassland and woodland; and woodland and forest. The inconsistency in the differences of taxa distribution in each vegetation type relative to the whole species pool might be explained by the number of species in our sample. It is important to note that accounting for geographical scale and species abundance are important to draw strong conclusions from phylogenetic community structure (Vamosi et al., 2009). The differences in forest and woodland phylogenetic structure suggest the species characterizing these environments evolved independently and based on Vorontsova et al. (2016), it appears that both grassland and woodland communities have a long evolutionary history. Hackel et al. (2018) demonstrated that the "Madagascar shade clade" diversified in Madagascar and nearby islands since the Miocene and that endemic C<sup>4</sup> grasses are the result of repeated dispersals to Madagascar mainly since the Late Miocene. It is likely that immigrating C<sup>4</sup> grasses already featured traits similar to those we observed, but further study and comparison with their non-Malagasy relatives is needed.

Our assessment of the grass community composition and their characteristics in grassland, woodland and forest shows that the woodland formation dominated by Uapaca bojeri is more likely a type of savanna and has been previously wrongly attributed as forest (Perrier de la Bâthie, 1921; Humbert, 1949; Koechlin et al., 1974; Gade, 1996). Our results support previous studies which distinguished forest and savanna based on tree functional traits and diversity (Ratnam et al., 2011; Dantas et al., 2013; Charles-Dominique et al., 2015). Uapaca bojeri, the dominant tree species in woodland does not occur in forest as defined here and our comparative analysis of bark thickness reinforces this conclusion. Uapaca bojeri (tapia) has bark at least twice the thickness of nearby forest species. Bark thickness is often presented as a key mechanism for tree survival in frequently burnt ecosystems (Gignoux et al., 1997; Hoffmann et al., 2003; Lawes et al., 2013; Pausas, 2015) where thick bark reduces the probability of stem death and insulates buds buried in stems (Hoffmann et al., 2003; Pausas, 2015). In addition, Uapaca bojeri has other traits characteristic of pyrophytic species (**Table 1**), including fire retardant leaves and epicormic resprouting (Kuhnholtz-Lordat, 1938). Our results support Kull's (2002) assertion that Uapaca bojeri dominated woodlands are shaped and maintained by fire, suggesting this woodland is a savanna rather than degraded forest. However, comparison of multiple tree traits related to fire resistance, tolerance and recruitment is needed. Previous research has found differences in specific leaf area and leaf:stem allocation patterns in tree species across savanna and forest boundaries (Hoffmann and Franco, 2003) and such trait comparisons along with assessments of conditions of plant recruitment could help further elucidate the functional ecology of these woodlands linked to the functional ecology of the ground layer. These woodlands, with numerous endemic species in both the over- and under-story, are part of the cultural landscape of Madagascar and have therefore undergone changes linked to wood collection and human use (Kull, 2002). Developing effective land management and conservation strategies of the unique and endemic flora of Madagascar requires development of an understanding of the functional characteristics of the species and vegetation types. Our study provides the first assessment of grass species composition and functional traits as being distinct between forest and the mosaic of grassland and Uapaca bojeri dominated woodlands, suggesting that functional characteristics of the ground layer can help delimit functionally distinct ecosystems. Further, the grasses of both grassland and woodland are phylogenetically highly diverse, consistent with our hypothesis that the evolutionary history of these ecosystems is certainly far older that human settlement which occurred only 6,000 years ago on Madagascar such that these woodlands of the central highlands are an ancient feature of these landscapes.

## AUTHOR CONTRIBUTIONS

CS, MV, and CL designed research. CS, MV, JH, GB, JW, and CL collected data. CS and JH analyzed data. CS, MV, and CL wrote the paper. JH, GB, SC, JW and VJ contributed to the interpretation and the revision of the work.

## FUNDING

CS was supported by KMCC and CERL by KMCC and University of Edinburgh. JH and GB received support from the University of Toulouse III through an Idex UNITI mobility grant (RECAC15EX110; UPS) and a DREIC project (VSR28AFRIQ; UPS). They are also supported by the French excellence project Labex TULIP (ANR-10-LABX-0041).

## REFERENCES


## ACKNOWLEDGMENTS

We thank T. Ramananantoandro and R. N. Razafinarivo from the forestry department of Antananarivo University for providing wood samples; British Ecological Society Travel Grant to Ecology Across Borders 2017 which enabled CLS to attend and present a poster at the conference; Madagascar National Parks (MNP), Direction générale des forêts (DGF) and Parc Botanique et Zoologique de Tsimbazaza (PBZT) for granting research permits, Dr. Tendro Tondrasoa (University of Aberdeen), Jan Kim (RBG Kew), and Tim Wilkinson (RBG Kew) for helping with data analysis, all Kew Madagascar Conservation Centre (KMCC) individual staff for support with permits and field trip, and all local community managers we worked with.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2018.00184/full#supplementary-material


FAO (2006). Global Forest Resources Assessment 2005. Rome: FAO.


and ecology. Bioinformatics 26, 1463–1464. doi: 10.1093/bioinformatics/ btq166


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Solofondranohatra, Vorontsova, Hackel, Besnard, Cable, Williams, Jeannoda and Lehmann. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Functional Traits of Trees From Dry Deciduous "Forests" of Southern India Suggest Seasonal Drought and Fire Are Important Drivers

Jayashree Ratnam<sup>1</sup> \*, S. K. Chengappa<sup>1</sup> , Siddarth J. Machado<sup>1</sup> , Nandita Nataraj <sup>1</sup> , Anand M. Osuri 2,3 and Mahesh Sankaran1,4

*<sup>1</sup> National Centre for Biological Sciences, Tata Institute of Fundamental Research, Bangalore, India, <sup>2</sup> The Earth Institute, Columbia University, New York, NY, United States, <sup>3</sup> The Nature Conservancy, Arlington, VA, United States, <sup>4</sup> School of Biology, University of Leeds, Leeds, United Kingdom*

#### Edited by:

*Daniel M. Griffith, Oregon State University, United States*

#### Reviewed by:

*Juan Ernesto Guevara Andino, Field Museum of Natural History, United States Kathleen M. Quigley, Michigan State University, United States*

> \*Correspondence: *Jayashree Ratnam jratnam@ncbs.res.in*

#### Specialty section:

*This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution*

> Received: *01 June 2018* Accepted: *09 January 2019* Published: *29 January 2019*

#### Citation:

*Ratnam J, Chengappa SK, Machado SJ, Nataraj N, Osuri AM and Sankaran M (2019) Functional Traits of Trees From Dry Deciduous "Forests" of Southern India Suggest Seasonal Drought and Fire Are Important Drivers. Front. Ecol. Evol. 7:8. doi: 10.3389/fevo.2019.00008* Two dominant biomes that occur across the southern Indian peninsula are dry deciduous "forests" and evergreen forests, with the former occurring in drier regions and the latter in wetter regions, sometimes in close proximity to each other. Here we compare stem and leaf traits of trees from multiple sites across these biomes to show that dry deciduous "forest" species have, on average, lower height: diameter ratios, lower specific leaf areas, higher wood densities and higher relative bark thickness, than evergreen forest species. These traits are diagnostic of these dry deciduous "forests" as open, well-lit, drought-, and fire-prone habitats where trees are conservative in their growth strategies and invest heavily in protective bark tissue. These tree traits together with the occurrence of a C<sup>4</sup> grass-dominated understory, diverse mammalian grazers, and frequent fires indicate that large tracts of dry deciduous "forests" of southern India are more accurately classified as mesic deciduous "savannas."

Keywords: functional traits, fire, savnnas, forests, southern India

## INTRODUCTION

Ecologists have long organized earth's vegetation into biomes or vegetation types that occur predictably under certain combinations of precipitation and temperature (Schimper, 1903; Holdridge, 1947; Walter, 1973; Whittaker, 1975). However, Whittaker (1975); recognized that there were certain climatic zones where vegetation was unpredictable; this zone of "ecosystems uncertain" includes large parts of the tropics where vegetation can be tropical forest, savanna, shrubland, or grassland, depending on seasonality, drought, and other disturbances (Whittaker, 1975; Bond et al., 2005; Bond, 2013). It is also now widely agreed that while temperature and precipitation together define the major biomes at continental scales, local scale turnover in vegetation types can be driven by sharp differences in underlying topographic and edaphic conditions (Ratter, 1992; Esler et al., 2015; Moncrieff et al., 2016; Miatto and Batalha, 2017) but more often by vegetation-disturbance feedbacks that can result in very different ecosystems within comparable climatic conditions (Bond et al., 2005; Staver et al., 2011; Hoffmann et al., 2012; Dantas Vde et al., 2013, 2016; Charles-Dominique et al., 2015; Pausas and Dantas, 2017).

Across the mesic tropics, where both wooded savannas and forests occur, sometimes in close proximity within a landscape, recent evidence clearly establishes a critical role for vegetation-fire feedbacks in maintaining these savannaforest transitions (Hoffmann et al., 2012; Dantas Vde et al., 2013; Charles-Dominique et al., 2015; Gray and Bond, 2015). Specifically, open lighted environments and frequent grassfuelled fires are associated with the savanna state, while closed, shaded environments and rare fires are associated with the forest state (Bond and Parr, 2010; Hoffmann et al., 2012). Reflecting these differences, the woody species that are characteristic of these different ecosystems are expected to vary in their morphological and physiological traits in ways that are adapted to their distinctive disturbance-environment regimes (Ratnam et al., 2011; Hoffmann et al., 2012; Dantas Vde et al., 2013).

Forest trees, which rarely encounter fire and are under selection to out-shade competitors are predicted to grow taller and develop wider canopies for a given diameter, and invest little in protective bark tissue (Ratnam et al., 2011). They may also have higher specific leaf area and lower wood density, which are supportive of rapid nutrient acquisition and growth (Poorter and Bongers, 2006; Chave et al., 2009; Miatto and Batalha, 2017). In contrast, savanna trees that grow in fire-prone environments are expected to invest heavily in thick bark that protects them from fire (Gignoux et al., 1997; Hoffmann et al., 2012; Dantas Vde et al., 2013). They may store much of their biomass in belowground tissue that is safe from fire, and aboveground, they may be shorter and have narrower canopies for a given diameter (Ratnam et al., 2011). They may also be expected to have lower specific leaf area and higher wood densities that support the more conservative growth strategies that are favored in their open and desiccating environments (Miatto and Batalha, 2017).

Here, we consider the functional traits of woody trees across different vegetation types in the southern Indian peninsula. Much of this region falls within the zone of "ecosystems uncertain" in the seasonal tropics (Bond, 2013, Figure 16.1); it supports a diversity of vegetation formations ranging from open thorn scrub to wooded grasslands to closed forests. However, a historical problem with the vegetation nomenclature in this region is that all the vegetation formations with some degree of woody cover are classified as "forests" (see Champion and Seth, 1968; Ratnam et al., 2016). Thus, both open and closed woody formations fall within this nomenclature such that savannas are not distinguished from forests. In past work, we have argued that large tracts of the vegetation type in this region classified as "dry deciduous forest" that are characterized by relatively open tree canopies in grassy understories are more correctly viewed as mesic savannas, while vegetation types with dense tree cover and non-grassy understories are true forest formations (Ratnam et al., 2011, 2016; Sankaran and Ratnam, 2013). Here, we extend these arguments to an examination of the traits of tree species in these habitats. If, as we argue, many of the dry deciduous "forests" in southern India are in fact deciduous "savannas" that are adapted to seasonal water-stress and frequent grass-fuelled fires (Ratnam et al., 2016), which is not the case for true forests, we expect traits of tree species characteristic of these habitats to reflect this difference.

## METHODS

We sampled dry deciduous and evergreen forests across seven sites in southern India (see inset in **Figure 1**). Rainfall across our sampling points in dry deciduous forests ranged from 516 to 1,260 mm, while rainfall across sampling points in evergreen forests ranged from 1,078 to 5,660 mm. Rainfall across this entire region is monsoonal such that all sites experience long dry seasons of 5–8 months per year (**Figure 1**).

We collected functional trait data from 1,350 individual trees of 75 dry deciduous "forest" (DD) species and 92 (EG) evergreen forest species (**Supplementary Table 1**), with 9 species common to both habitats. Specifically, we measured specific leaf area (SLA), tree height, girth at breast height (GBH), wood density, and bark thickness. For each species, in each habitat, we collected traits from multiple individuals (3–25, depending on the commonness of the species).

For SLA, five mature sun leaves were collected from the canopy of each individual and scanned in a flat-bed scanner on the same day. Leaf areas of scanned leaves were calculated using either Blackspot (Varma and Osuri, 2013) or Image J software (Rasband, 2014). Leaves were then oven-dried and dry weights measured in a precision balance. SLA was calculated as leaf area per unit dry weight. GBH was measured with a tape-measure, 1.3 m above the ground. Tree heights were measured using a clinometer or, alternately, a laser rangefinder when the observer could stand at the base of the trunk and sight the top of the tree with the laser beam. Wood cores were collected using an increment borer, and placed into sealed plastic bags with moist cotton. Back in the field station, they rehydrated in water for 1 h. Fresh wood volume was estimated by water displacement, following which samples were oven-dried at 65◦C for 72 h and weighed, and wood density estimated as oven-dried weight divided by fresh volume (Chave, 2005). Bark thickness was measured at breast height on the trunk, using a bark gauge. Using the sharp end of the gauge, we gently gouged through the outer layers of bark until it gave way to the thick inner wood which was visibly different tissue, and measured bark thickness to this point. Relative bark thickness was calculated as: (Bark thickness<sup>∗</sup> 2/Diameter)<sup>∗</sup> 100.

We used Generalized Linear Mixed Models (GLMM) to look at how traits differed between the two vegetation types. Vegetation type (DD, EG) was included as a fixed effect, with species nested within sites as the random effect. Analyses were carried out using the "lme4" package, and the "car" package was used to conduct Type II Wald Chi square tests to assess statistical significance of the fixed effects (Bates et al., 2017). Values for SLA, H:D, and relative bark thickness were log transformed to meet model assumptions.

## RESULTS

Consistent with our expectations, dry deciduous "forest" species differed significantly from evergreen forest species for the four traits measured. DD species had, on average, lower

height: diameter ratios (Mean ± SE, DD: 0.360 ± 0.01, EG: 0.576 ± 0.017, χ <sup>2</sup> = 184.71, df = 1, p < 0.0001) higher relative bark thickness (DD:8.75 ± 0.47%, EG:4.26 ± 0.26%, χ <sup>2</sup> = 50.03, df = 1, p < 0.0001), lower SLA (DD: 94.58 ± 2.42 g/cm<sup>2</sup> , EG:106.31 ± 2.96 g/cm<sup>2</sup> , χ 2 = 5.801, df = 1, p = 0.016), and marginally higher wood densities (DD: 0.619 ± 0.015 g/cm<sup>3</sup> , EG:0.603 ± 0.011 g/cm<sup>3</sup> , χ <sup>2</sup> = 5.87, df = 1, p = 0.015) than EG species (**Figure 2**).

## DISCUSSION

Our cross-site comparison of the functional traits of tree species from dry deciduous and evergreen forests of southern India indicate that dry deciduous "forests" are fundamentally different in the environments and species traits that they support than are evergreen forests, even though they occur in close proximity in many regions. Dry deciduous "forest" trees have lower height to diameter ratios, produce leaves with lower specific leaf areas, and

have higher wood densities. Taken together, this suite of traits is suggestive of slower and more conservative growth strategies in these trees than in their evergreen forest counterparts (Poorter and Bongers, 2006; Chave et al., 2009; Miatto and Batalha, 2017). Such strategies may indicate either relatively lower resource environments, or less competitive ones (Miatto and Batalha, 2017). Clearly, dry deciduous "forests" occur in areas with less rainfall than evergreen forests, although the two forest types overlap at the wetter end of the dry deciduous zone. Seasonal water stress is thus likely to be more pronounced in these habitats. With relatively lower tree densities, dry deciduous "forests" may also be both more desiccating, and less competitive environments for adult trees that grow in them.

Critically, dry deciduous "forest" trees have dramatically higher investment in bark, producing, on average, almost twice as much bark for a given diameter as evergreen forest trees. Thick barks are known to serve multiple protective functions including drought resistance, fire resistance, and prevention of stem damage (Poorter et al., 2014; Schafer et al., 2015). For fire-prone environments, there is overwhelming evidence for a primary role of thick barks in enabling trees to persist through fires (Lawes et al., 2011; Hoffmann et al., 2012; Pausas, 2015; Pellegrini et al., 2017) and several recent analyses across the globe have shown that fire-resistant savanna trees are characterized by thicker bark than fire-sensitive forest trees (Hoffmann et al., 2012; Dantas Vde et al., 2013; Lawes et al., 2013; Pellegrini et al., 2017). Thicker barks may also increase drought resistance in drier environments (Poorter et al., 2014; Schafer et al., 2015) as in these dry deciduous "forests," but this function may be secondary to a strategy of deciduousness where species avoid drought by being leafless in the dry season.

Several analyses of contemporary fire regimes in the subcontinent confirm that dry deciduous "forests" burn much more frequently (fire return intervals range from 1 to 6 years) and extensively (10–50% of these landscapes burn annually) than do other forest types (Kodandapani et al., 2008; Kodandapani, 2013; Srivastava and Garg, 2013; Mondal and Sukumar, 2016; Reddy et al., 2017). Together with the trait comparisons we report in this study, which are consistent with predicted trait differences between mesic savanna and forest species (Ratnam et al., 2011), these data suggest that large tracts of dry deciduous "forests" in this region, characterized by trees in grass-dominated understories, are in fact deciduous "savannas," where seasonal water-stress and fire are important drivers. Further, these habitats are also home to an ancient and diverse mammalian herbivore assemblage, ranging from smaller-bodied spotted deer to largebodied mega-herbivores like the Asian bison and elephant (Bibi and Métais, 2016; Sankaran and Ahrestani, 2016). Large mammalian grazers, which are supported by the forage in grassdominated understories are characteristic of savannas globally, additional evidence that these "forests" are indeed "savannas."

These confirmatory findings have important implications for how fire is managed in these ecosystems. At the current time, these habitats are managed under a blanket policy of fire exclusion and suppression and fire-setting is a punishable offense according to the Indian Forest Act (Ratnam et al., 2016; Thekaekara et al., 2017). These policies, which stem from a historical misreading of these ecosystems as "forests" are inimical to the conservation and sustainable management of these savannas (Lehmann and Parr, 2016; Griffith et al., 2017). First, the notion that fires are always undesirable disturbances in these ecosystems, widespread amongst both managers and vegetation scholars, prevents nuanced understanding, and appropriate research on the ecology of these systems. Second, these policies preclude the use of fire as a management tool, an opportunity lost in many areas where fire exclusion has been associated with severe non-native shrub invasions (Sundaram et al., 2012) and associated losses of native biodiversity (Ramaswami and Sukumar, 2011; Sundaram and Hiremath, 2012). Official recognition of these tree-grass ecosystems as drought- and firedriven mesic savannas will provide multiple opportunities for vital research and effective management.

## AUTHOR CONTRIBUTIONS

JR led and wrote the paper. MS and JR planned the trait comparisons. SKC, SJM, NN, AMO, JR, and MS collected data. All authors commented on the manuscript.

## FUNDING

Funding for this study came from Department of Science and Technology, Government of India (SERB/SR/SO/PS/78/2012), core funds from NCBS-TIFR to MS and Rufford Small Grants Foundation (11086-2) to AMO.

## ACKNOWLEDGMENTS

We thank the Forest Departments of Karnataka, Tamil Nadu, Telangana, and Andhra Pradesh for permits to

## REFERENCES


conduct this study. We thank Dayani Chakaravarthy, Chintan Seth, VijayKumar S, Swapna Nellaballi, and Catherin Arockia for help with data collection and processing.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2019.00008/full#supplementary-material


Rasband, W. S. (2014). ImageJ. Bethesda, MD: U. S. National Institutes of Health.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Ratnam, Chengappa, Machado, Nataraj, Osuri and Sankaran. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Inserting Tropical Dry Forests Into the Discussion on Biome Transitions in the Tropics

Kyle G. Dexter 1,2 \*, R. Toby Pennington2,3, Ary T. Oliveira-Filho<sup>4</sup> , Marcelo L. Bueno<sup>5</sup> , Pedro L. Silva de Miranda<sup>1</sup> and Danilo M. Neves 4,6

<sup>1</sup> School of GeoSciences, University of Edinburgh, Edinburgh, United Kingdom, <sup>2</sup> Royal Botanic Garden Edinburgh, Edinburgh, United Kingdom, <sup>3</sup> Geography, University of Exeter, Exeter, United Kingdom, <sup>4</sup> Departamento de Botânica, Universidade Federal de Minas Gerais, Belo Horizonte, Brazil, <sup>5</sup> Laboratório de Ecologia e Evolução de Plantas, Universidade Federal de Viçosa, Viçosa, Brazil, <sup>6</sup> Department of Ecology and Evolutionary Biology, University of Arizona, Tucson, AZ, United States

Tropical moist forests and savannas are iconic biomes. There is, however, a third principal biome in the lowland tropics that is less well known: tropical dry forest. Discussions on responses of vegetation in the tropics to climate and land-use change often focus on shifts between forests and savannas, but ignore dry forests. Tropical dry forests are distinct from moist forests in their seasonal drought stress and consequent deciduousness and differ from savannas in rarely experiencing fire. These factors lead tropical dry forests to have unique ecosystem function. Here, we discuss the underlying environmental drivers of transitions among tropical dry forests, moist forests and savannas, and demonstrate how incorporating tropical dry forests into our understanding of tropical biome transitions is critical to understanding the future of tropical vegetation under global environmental change.

#### Edited by:

Colin Osborne, University of Sheffield, United Kingdom

### Reviewed by:

Imma Oliveras, University of Oxford, United Kingdom Liam Jude Langan, Senckenberg Biodiversität und Klima Forschungszentrum (SBiK-F), Germany

\*Correspondence:

Kyle G. Dexter kyle.dexter@ed.ac.uk

#### Specialty section:

This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution

> Received: 01 May 2018 Accepted: 29 June 2018 Published: 24 July 2018

#### Citation:

Dexter KG, Pennington RT, Oliveira-Filho AT, Bueno ML, Silva de Miranda PL and Neves DM (2018) Inserting Tropical Dry Forests Into the Discussion on Biome Transitions in the Tropics. Front. Ecol. Evol. 6:104. doi: 10.3389/fevo.2018.00104 Keywords: tropical dry forest, tropical moist forest, savanna, biomes, fire, soil fertility, water stress, deciduousness

## INTRODUCTION

Predicting vegetation change in the tropics depends on understanding the drivers of transitions among major vegetation types, or biomes. Climatic factors, such as mean annual precipitation (MAP) and its seasonality are of obvious importance, but edaphic factors, and disturbance via fire, humans and herbivores also play key roles. Recent large-scale studies across the tropics have focused on transitions between forest and savanna (Hirota et al., 2011; Staver et al., 2011; Oliveras and Malhi, 2016; Xu et al., 2016; Langan et al., 2017). While there is value to simplifying vegetation concepts in the tropics, we believe the simplification used by these authors in defining "forest" goes perhaps one step too far. There are two principal kinds of forest in the lowland tropics, moist forests and dry forests. With very few exceptions (e.g., Hirota et al., 2010; Lehmann et al., 2011), studies of biome transitions in the tropics have either failed to distinguish them, or have completely ignored dry forests, focusing solely on moist forests when using the term "forest." The aim of this perspective is to discuss biome transitions in the tropics and their underlying drivers, while including dry forests in the discussion. We focus on transitions among savanna, moist forest and dry forest, the three biomes in the lowland tropics with a substantial tree component.

Tropical moist forest and savanna are relatively well understood at a global scale compared to tropical dry forests (Pennington et al., 2018). Moist forest is tall, multi-stratal, and with a closed canopy. Tropical moist forests include tropical rain forests, as well as forests with lower rainfall where soil moisture is maintained throughout the year, via edaphic factors such as proximity to rivers, or water recycling, allowing most trees to be evergreen (Guan et al., 2015). The understory is often dominated by saplings of taller-statured tree species, although small tree and shrub species are present. Terrestrial forbs and grasses are a minor component of diversity and biomass. Savanna is a more open environment, where tree species are present, but individuals do not form a closed canopy. There is a significant understory grass component, which is flammable, and fires are common. Tree species that occur in savannas are adapted to these recurring fires (Simon and Pennington, 2012), and regular fires are necessary for the maintenance of savanna biodiversity (Parr et al., 2014; Durigan and Ratter, 2016; Abreu et al., 2017). Some savannas in the paleotropics (e.g., miombo woodlands in Africa, deciduous dipterocarp forests in southeast Asia) are generally referred to as dry forests, but we consider them as savannas given that they have a grassy understory and experience regular fire (Ratnam et al., 2011; Dexter et al., 2015; Pennington et al., 2018).

Tropical dry forests vary greatly in structure, from tall, closed canopy forest to short scrub vegetation, occasionally not forming a closed canopy, especially in drier areas (Pennington et al., 2000). They are distinct from savanna in not having a significant grass component and not experiencing regular fires (Murphy and Lugo, 1986; Gentry, 1995). In fact, regular fires would be lethal for many of the characteristic life forms and taxa of tropical dry forest (e.g., cacti; Mooney et al., 1995). This is not to say that dry forests never experience fires. Even moist forests can experience fire under extreme drought conditions (Aragão et al., 2016). Rather, damaging fire is sufficiently rare in dry forests such that fire-intolerant species can persist in the landscape as metapopulations (Hanski, 1998). The exact threshold of fire return interval or intensity involved in the tropical dry forest savanna transition is poorly understood and likely to vary with the broader environmental context (e.g., soil fertility and annual precipitation; Hoffmann et al., 2012; Murphy and Bowman, 2012), and we suggest that this should be a priority for further study. Tropical dry forest is distinct from moist forest in its seasonal drought stress, which leads many tree species to lose their leaves in the dry season (Reich and Borchert, 1984; Murphy and Lugo, 1986). The combination of seasonal drought stress and lack of fire leads to ecosystem function in dry forests that is markedly different from savannas or moist forests, which justifies their distinction as a unique biome.

## BIOMES IN LOWLAND TROPICAL SOUTH AMERICA

We focus this review on continental lowland tropical South America (LTSA), where we have conducted most of our research. In a recent study (Silva de Miranda et al., in press), we used an unsupervised classification, or hierarchical clustering, of sites based on their tree species composition (see inset in **Figure 1**), followed by interpretation of the resulting cluster using site information on vegetation physiognomy (savanna vs. forest) and leaf flush regime (evergreen vs. semideciduous vs. deciduous) to delimit and map biomes across LTSA east of the Andes (**Figure 1**). Moist forests fell in two major groups in the cluster and occurred in two large geographic blocks, one in the Amazon basin and another along the Atlantic coast of Brazil. Semideciduous forests did not form a distinct group in the cluster and were instead mixed with evergreen, moist forest sites. Semideciduous forests are often found in drier regions, where they occur along rivers, lake margins and submontane areas with orographic precipitation. Savanna formed a single group in the cluster and was most prevalent in central Brazil, in an area commonly termed the Cerrado.

Tropical dry forest also formed a single group in the cluster, which was comprised almost entirely of forest sites with deciduous phenology. In LTSA, the largest block of dry forest occurs in the Caatinga region of northeast Brazil. The Caatinga has been referred to as a biome (Hirota et al., 2010; Instituto Brasileiro de Geografia e Estatística, 2012), although as a region it contains non-dry forest habitat (e.g., patches of savanna). Further, tropical dry forest is found outside of the Caatinga, in patches throughout the Cerrado, in an area spanning the Pantanal and Chiquitania in Bolivia (**Figure 1**), and scattered more widely across the Neotropics (DRYFLOR, 2016). While Silva de Miranda et al. (in press) broadly assessed climatic overlaps amongst biomes, they did not focus on the environmental drivers of transitions between individual biomes. That is the goal of the present manuscript.

## TRANSITIONS BETWEEN TROPICAL SAVANNA AND DRY FORESTS

A common view of dry forests in the tropics is that they are transitional between savannas and moist forests along precipitation gradients (e.g., Whittaker, 1970; Malhi et al., 2009). If, however, we examine how the sites featured in **Figure 1** are distributed over variation in MAP (**Figure 2**), a more complex picture emerges. Moist forests do occur under wetter conditions than savanna, but dry forests are largely found under drier conditions. Below 1,000 mm MAP, savanna quickly disappears and dry forest becomes the only tree-dominated vegetation type. The largest area of these arid dry forests is found in the Caatinga region of northeast Brazil. As discussed above, a key distinction between savanna and dry forest is the regularity of fire, and in these dry conditions, there is not sufficient biomass build-up to sustain regular fires (Van Der Werf et al., 2008). In particular, this reflects the relative lack of grasses. The tree species that are present are able to tolerate severe drought, but they do not invest in adaptations for fire, such as thick bark or underground stems, characteristic of savanna species (Simon and Pennington, 2012).

There is also extensive occurrence of the dry forest biome under the same precipitation conditions as savanna. These are dry forests found in the Cerrado region and around the Pantanal and Chiquitania regions of Brazil and Bolivia. Within the Cerrado, dry forests are known to occur on and around calcareous outcrops, where soils have higher phosphorus and base cation concentrations (Ratter et al., 1978; Furley and Ratter, 1988; Oliveira-Filho and Ratter, 2002; Neves et al., 2015). On these soils, trees can grow more quickly, have better chances of escaping the "fire trap" and are more likely to

form a closed canopy (Hoffmann et al., 2009). There can be positive feedbacks between tree growth, grass exclusion and fire mitigation that leads to a forested vegetation (Hoffmann et al., 2012; Silva et al., 2013; Pausas and Dantas, 2017). Calcareous outcrops in the Cerrado also have poorly developed, shallow soils, and vegetation occurring on them may experience greater drought stress than surrounding vegetation, thus making them similar to the arid dry forests, which lack fire because of insufficient fuel build-up. However, it is likely that soil fertility is relevant for the presence of dry forests around calcareous outcrops as a different vegetation, cerrado rupestre, which is floristically related to savanna vegetation, is found on noncalcareous outcrops in the Cerrado (Ribeiro and Walter, 1998). Whichever factor is more important (soils with high fertility or low water-holding capacity), it is evident that the same

et al. (in press), but are distinguished here. We exclude sites south of 23◦S latitude or above 1,000 m elevation.

drought-tolerant, fire-intolerant tree species and lineages that dominate vegetation in the arid Caatinga are also found in dry forest patches in the moister Cerrado (Prado and Gibbs, 1993; Neves et al., 2015; DRYFLOR, 2016; Silva de Miranda et al., in press).

Dry forest and savanna vegetation also intermingle in the Chiquitania region of Bolivia, but here dry forest predominates and savanna occurs in patches, which may be because soils in the Chiquitania are more fertile on average than in the Cerrado (Silva de Miranda et al., in press). Some of the savannas that are present in the Chiquitania region may represent dry forest that has been degraded by logging or anthropogenic fire (Devisscher et al., 2016), highlighting that human land-use patterns can readily drive transitions between dry forest and savanna. However, "oldgrowth savannas" that would exist independent of anthropogenic

influence are also clearly present in Bolivia (Power et al., 2016; Veldman, 2016).

If fire is excluded from savanna vegetation, it first converts to a formation with a higher percentage of tree canopy cover termed "cerradão" (Durigan and Ratter, 2006, 2016), which translates from the Portuguese as "big cerrado" and is generally considered as a forest. Cerradão shares some tree species with tropical dry forest and comparatively few with semideciduous and evergreen moist forests (Bueno et al., 2018), even though the latter are present near savanna vegetation along river courses and lake margins that have year-round water availability (Ribeiro and Walter, 1998). Tree species from semideciduous and evergreen forests may be less likely to immigrate into cerradão than typical dry forest tree species because they are not adapted to seasonal drought. While cerradão is initially comprised of fire-adapted tree species from the cerrado, dry forest tree species colonize this environment if propagules are available, fire remains absent and soils are sufficiently fertile. These dry forest tree species may eventually outcompete cerrado tree species, since they do not invest in fire defense (Ratajczak et al., 2017), and in the prolonged absence of fire, cerradão may transition to a dry forest if there are positive feedback cycles between forest vegetation, lack of fire and soil fertility (Silva et al., 2013; Pellegrini et al., 2014, 2018). However, if the underlying soils remain poor and/or if there are high aluminum concentrations in the soil that do not attenuate over time, then cerrado tree species, which are adapted to infertile, aluminum-rich soils may continue to dominate the vegetation.

Overall, the savanna-dry forest transition is distinct from the savanna-moist forest transition in two key ways: (1) the contrasting role of water availability (lowest in dry forest, intermediate in savanna and highest in moist forest) and (2) the potentially critical importance of soil fertility for the savannadry forest transition (savanna and moist forests are similar in generally having infertile, acidic soils).

## TRANSITIONS BETWEEN MOIST FOREST AND DRY FOREST

Both of these biomes are forest, but they function in distinct ways. Dry forests are found in areas with marked precipitation seasonality, which leads most species to lose their leaves during the dry season and has significant implications for nutrient cycling (Reich and Borchert, 1984; Murphy and Lugo, 1986). Many moist forests also experience seasonality in water availability (e.g., in the southern and eastern Amazon), but the dry season is three months or less and subsurface water remains available to trees (Guan et al., 2015). The systems also differ in the rate at which they accrue and cycle carbon, with trees in moist forests growing more quickly and storing more total carbon (Murphy and Lugo, 1986; Poorter et al., 2017). There are often differences in soil fertility, with dry forests occurring on more fertile soils, which facilitates their ability to shed their leaves as they can readily afford to grow new ones. However, high rainfall in moist forests results in nutrient leaching, and this correlation between soil fertility and biome identity may be due to overriding climatic factors (Webb, 1968; Hall and Swaine, 1976).

In LTSA, there are multiple areas of contact between the moist and dry forest biomes (**Figure 1**). One transition zone is in northeastern Brazil, where evergreen Atlantic forest on the coast transitions to dry forest in the arid Caatinga. In between the two lies a band of semideciduous forests. Another transition zone is found in the Chiquitania region of eastern Bolivia and adjacent areas of Brazil, where there is a gradual transition over 200+ km of geographic distance, largely covered by semideciduous forest. We suggest that transitions between dry forest and moist forest are primarily mediated by water availability and that intermediate states are possible in zones of intermediate water availability (Oliveira-Filho and Fontes, 2000; Oliveira-Filho et al., 2006). This contrasts with transitions between savanna and moist forest that can be more abrupt and may represent alternative stable states (although see Lloyd and Veenendaal, 2016).

## SEMIDECIDUOUS FORESTS

Previous studies have variously grouped semideciduous forests with the dry forest biome (Murphy and Lugo, 1986; Pennington et al., 2000) or the moist forest biome (DRYFLOR, 2016; Silva de Miranda et al., in press). In fact, as discussed above, these forests may be transitional between the two. Semideciduous forests in LTSA have few endemic tree species and instead contain tree species associated with the dry or moist forest biomes (Oliveira-Filho and Fontes, 2000). As moist forests contain many more tree species than dry forests (Esquivel-Muelbert et al., 2017), we suggest that they may contribute more species to semideciduous forests simply via mass effects (Shmida and Wilson, 1985), and this may be why they group with moist forests in clustering analyses based on presence versus absence of tree species (as in Silva de Miranda et al., in press). If abundance information were to be taken into account (e.g., via inventory plot data), we hypothesize that semideciduous forests may cluster with moist or dry forest based on the proportion of individuals belonging to moist versus dry forest tree species. It is clear that future comparative studies across dry, moist and semideciduous forests are needed to understand their origins and how they compare in terms of ecosystem function. Their geographically variable species composition and lack of endemic species suggests that semideciduous forests may have been independently and recently assembled in different ecotonal areas.

## BIOME TRANSITIONS TO DRY FOREST OUTSIDE THE NEOTROPICS

As in South America, moist forest—savanna transitions have been studied extensively on other continents, but transitions to dry forest have received less attention. This is partly because it is unclear where dry forest exists outside of the Neotropics (Lock, 2006; Dexter et al., 2015; Pennington et al., 2018). In a recent study, Linder (2014) delimited and described the main "floras" of Africa, which are large-scale units of vegetation that have a distinct evolutionary and biogeographic history and differ in their present-day plant taxonomic composition. Linder did not assign the term "biome" to these vegetation units, although his "floras" correspond to several previously defined biomes. There is a "savanna flora," which readily corresponds to the savanna biome, and a "lowland forest flora" that largely corresponds to the moist forest biome. Linder postulated an "arid flora" that is most evident in the Horn of Africa (Somalia, Ethiopia and northern Kenya), but is also present in arid regions of Angola and Namibia. The distribution of this "arid flora" largely overlaps the distribution of the "succulent biome" in Africa, as proposed by Schrire et al. (2005). The "arid flora" or "succulent biome" is similar to dry forest in arid regions of the Neotropics in that there is not adequate water availability to allow for sufficient biomass build-up to sustain regular fires. Thus, as in the Neotropics, water availability may be one environmental factor that underlies transitions between savanna and dry forest in Africa.

Soil fertility is another significant factor that has been shown to underlie savanna-dry forest transitions in the Neotropics. An important question for future research in Africa is whether, amongst its great expanses of savanna, there is distinct vegetation that does not regularly burn, is found on more fertile soils and shows greater floristic similarity with the "arid flora" of Linder (2014) than it does with surrounding savanna vegetation. It may be that in Africa, a higher abundance of large herbivores, including elephants, favors grasses over trees, leading to a more open savanna vegetation with more frequent fires, even in areas of higher soil fertility (Charles-Dominique et al., 2016; Pellegrini et al., 2017). If dry forests are not found on fertile soils in more mesic areas of Africa, there may not be moist forest-dry forest transitions on this continent, because the areas mapped as belonging to the "arid flora" or "succulent biome" are completely separated from moist forest regions by large areas of savanna (Schrire et al., 2005; Linder, 2014).

In the tropical regions of continental Asia and in Malesia, moist forest is the predominant vegetation type, although drier forest formations are present (e.g., deciduous forests in the Western Ghats and dry dipterocarp forests in Indochina). As we have discussed elsewhere (Dexter et al., 2015; Pennington et al., 2018), the majority of these drier forest formations have a significant grassy component in the understory, burn regularly and may be better considered as savannas (Ratnam et al., 2011). The succulent biome of Schrire et al. (2005) is mapped as present in arid regions of northwest India and extending across the coast of Pakistan and Iran to the Arabian Peninsula. Thus, an arid form of the dry forest biome may occur in Asia, as in Africa, and water availability may underlie transitions between vegetation in seasonally dry areas that regularly burns (what we term savanna) and that does not regularly burn (what we term dry forest). As with Africa, future research in Asia should assess if there are vegetation formations in seasonally dry, yet not arid, areas that: (1) are found on fertile soils, (2) do not regularly burn and (3) show greater floristic similarity with arid areas than with surrounding vegetation that does regularly burn. This will help determine if soil fertility is also important in understanding savanna-dry forest transitions in Asia, as it is in the Neotropics.

## CONCLUSIONS

The aim of this perspective has been to bring the tropical dry forest biome into discussions of biome transitions in the tropics. Previous studies of tropical biome transitions have largely focused on forest-savanna transitions, with all forests being considered as a single biome. In fact, there are many kinds of forests in the tropics, some of which are distinct from each other in species composition and ecosystem function and represent different biomes (i.e., dry vs. moist forest) and others which are more difficult to classify (e.g., semideciduous forests, cerradão). Water availability is a key factor underlying tropical biome transitions. While forests are often thought to occur under wetter conditions than savannas, tropical dry forest is actually more prevalent in areas of lower water availability (<1,000 mm MAP). Meanwhile, soil fertility, which has received limited attention in studies of biome transitions, is also critical in the Neotropics, and merits future research on other continents. More generally, recognizing tropical dry forest as a distinct biome within the tropics should improve the accuracy of modeling studies that aim to predict the future of tropical vegetation and ecosystem function under global environmental change.

## AUTHOR CONTRIBUTIONS

KD and RP wrote the first draft of the manuscript and contributed to revising the manuscript. AO-F, MB, PS and DN contributed to revising and improving the manuscript.

## REFERENCES


## FUNDING

KD, RP, and DN thank the Natural Environment Research Council (UK) for funding via NE/I028122/1. KD thanks the Leverhulme Trust for supporting him during the time this study was completed via an International Academic Fellowship. PS thanks the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior–Brazil (CAPES) for support via the Science without Borders Programme (grant 99999.013197/ 2013-04).


Savanna," in The Cerrados of Brazil, eds P. S. Oliveira and R. J. Marquis (New York, NY: Columbia University Press), 91–120.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Dexter, Pennington, Oliveira-Filho, Bueno, Silva de Miranda and Neves. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Transplant Experiments Point to Fire Regime as Limiting Savanna Tree Distribution

#### Nicola Stevens <sup>1</sup> \*, Sally A. Archibald<sup>2</sup> and William J. Bond3,4

*<sup>1</sup> Global Change and Biodiversity Group, Department of Botany and Zoology, Stellenbosch University, Stellenbosch, South Africa, <sup>2</sup> School of Animal, Plant and Environmental Sciences, University of the Witwatersrand, Johannesburg, South Africa, <sup>3</sup> Department of Biological Sciences, University of Cape Town, Cape Town, South Africa, <sup>4</sup> Fynbos Node, South African Environmental Observation Network, Cape Town, South Africa*

Plant species range shifts are predicted to occur in response to climate change. The predictions are often based on the assumption that climate is the primary factor limiting the distribution of species. However the distribution of grassy biomes in Africa cannot be predicted by climate alone, instead interactions between vegetation, climate and disturbance structure the ecosystems. To test if climatic variables, as predicted by an environmental niche model, determine the distribution limits of two common savanna tree species we established a transplant experiment at a range of latitudes and altitudes much broader than the distribution limits of our study species. We planted seedlings of two common savanna trees, *Senegalia nigrescens* and *Colophospermum mopane,* at eight paired high and low elevation sites across an 850 km latitudinal gradient in South African savannas. At each site seedlings were planted in both grassy and cleared plots. After 2 years of growth, rainfall, temperature and location inside or outside their distribution range did not explain species success. Grass competition was the only variable that significantly affected plant growth rates across all sites, but grass competition alone could not explain the distribution limit. Species distributions were best predicted when maximum tree growth rates were considered in relation to local fire return intervals. The probability of sapling escape from the fire trap was the most likely determinant of distribution limits of these two species. As trees grew and survived 100 s of kilometers south of their current range limits we conclude that climate alone does not explain the current distribution of these trees, and that climate change adaptation strategies for savanna environments based only on climatic envelope modeling will be inappropriate.

Keywords: species distribution, range limits, savanna, range shifts, transplant, grass competition, *Colophospermum mopane*, *Senegalia nigrescens*

## INTRODUCTION

Global change is transforming the planet (Ellis, 2011; Scheffers et al., 2016). One of the most pervasive impacts, which has severe economic and ecological impacts, is that of shifting plant and animal ranges (Parmesan and Yohe, 2003; Thomas et al., 2004; Chen et al., 2011). The shifts are most frequently attributed to the warming of the earth, with the expectation that plant and animal distributions will shift toward the cooler poles or higher elevations (Lenoir et al., 2008; Chen et al., 2011). These expectations (Hutchinson, 1957; Chase and Leibold, 2003) are often derived from the

#### *Edited by:*

*Daniel M. Griffith, Oregon State University, United States*

#### *Reviewed by:*

*Robert Bagchi, University of Connecticut, United States Katy Morgan, University of Bath, United Kingdom*

> *\*Correspondence: Nicola Stevens nicolastvns@gmail.com*

#### *Specialty section:*

*This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution*

*Received: 27 March 2018 Accepted: 23 August 2018 Published: 18 September 2018*

#### *Citation:*

*Stevens N, Archibald SA and Bond WJ (2018) Transplant Experiments Point to Fire Regime as Limiting Savanna Tree Distribution. Front. Ecol. Evol. 6:137. doi: 10.3389/fevo.2018.00137* Hutchinson concept of the fundamental and realized niche. The fundamental niche describes the set of environmental factors that determine the space where a species could occur in the absence of biotic interactions (Hutchinson, 1978; Chase and Leibold, 2003). The realized niche describes where a species actually occurs after accounting for the effects of inter-specific competition. This concept has been applied almost globally as a framework to describe the geographic distribution of species (Gaston, 2003; Thomas, 2010). However this framework is generally simplified to assume that climate shapes a plant's distribution range, and is the basis of attempts to predict the future ranges of plants using the correlative climate envelope approach (Thomas et al., 2004; Hijmans et al., 2005) where the current distribution of a species is mapped in climate space and any future shifts in climate space are expected to produce an accompanying shift in the future species distribution range (Elith and Leathwick, 2009). Although there is debate as to the broad applicability of climate models, (Davis et al., 1998; Araújo and Luoto, 2007; Svenning and Skov, 2007; Araújo and Peterson, 2012), the underlying general assumption that climate alone is the primary determinant of species distribution limits (Holdridge, 1967; Gaston, 2003) is seldom questioned.

Much of the empirical support for these climate based assumptions comes from studies in temperate and boreal northern hemisphere systems (Parmesan and Yohe, 2003; Thomas, 2010) where temperature is often the primary factor that limits plant production (Jolly et al., 2005; Körner and Basler, 2010).There remains an acute absence of detailed studies on distribution limits of plants in the tropics (Blach-Overgaard et al., 2010; Thomas, 2010). Here, authors acknowledge that the determinants of ranges in tropical regions are more complex (Root et al., 2003; Blach-Overgaard et al., 2010; Thomas, 2010). Water availability and biotic interactions are considered increasingly important in these systems (Woodward, 1987; Thomas et al., 2004; Jolly et al., 2005) but the underlying assumption that the distribution ranges can be predicted by climate alone remains common. However, we contend that this underlying principle and approach is not sufficient to predict the current, or future distribution of plants occurring in Africa's grassy biomes.

In tropical savanna and grassland biomes, interactions between fire, mammal herbivory and climate structure savannas (Bond et al., 2005; Staver et al., 2011; Lehmann et al., 2014). Tropical grassy biomes (TGBs) are therefore considered the most extensive vegetation area not at equilibrium with climate (Bond et al., 2005; Pausas and Ribeiro, 2017). Fire and mammal herbivory are dominant environmental filters that define community assemblages (Staver et al., 2012; Charles-Dominique et al., 2016; Pausas and Ribeiro, 2017), with the relative importance of these processes changing across the rainfall gradient (Archibald and Hempson, 2016). Herbivory can cause local species extinctions, alter competitive interactions, and reduce or enable tree recruitment (Prins and van der Jeugd, 1993; Bond and Loffell, 2001; Augustine and Mcnaughton, 2004; Asner et al., 2009), and fire regimes can alter the community structure and composition (Trapnell, 1959; Higgins et al., 2007; Pausas and Ribeiro, 2017). Whilst these processes can act on woody plants at all demographic stages (Midgley and Bond, 2001), the demographic stages which are most impacted are the seedling establishment and sapling stage where a critical population bottleneck occurs (Higgins et al., 2000). Here, regular disturbance traps saplings in the grass layer, where they are vulnerable to further topkill by fire and herbivory (Higgins et al., 2000; Midgley et al., 2010; Werner and Prior, 2013). Importantly these suppressed saplings are not reproductive and repeated topkill can prevent the recruitment of reproductive individuals, which would eventually lead to local extinction (Higgins et al., 2000). When trees exceed ∼3 m in height, they will often be sufficiently fire proof to be released from the fire trap (Bond, 2008; Wakeling et al., 2012). Therefore the ability of a plant to escape this trap to reach reproductive maturity within the prevailing climatic and disturbance context is the critical component that will determine the presence of a plant in a region.

To explore potential climatic constraints within a disturbance limited system we tested how well the traditional climate based framework of plant species distribution can explain patterns of species distribution ranges in consumer controlled systems. We established a transplant experiment across a latitudinal gradient across South African savannas inside and outside the distribution ranges of two common savanna tree species; Colophospermum mopane (Kirk ex J.Léonard) and Senegalia nigrescens [(Oliv.) P. J. H. Hurter]. The transplant experiment design was guided by the outputs of species distribution models and literature which describe the climatic determinants for the study species; Colophospermum mopane. Low winter temperatures (below a minimum July average of 8◦C), high rainfall (greater than 800 mm), and the presence of frost and increases in rainfall are commonly reported climatic limits to the occurrence of this species (Henning and White, 1974; Burke, 2006; Makhado et al., 2014; Scholtz et al., 2014; Stevens et al., 2014b).

To explore the role of these proposed limits the experiment was established across a climate gradient across a wide range of summer rainfall and temperature conditions, whilst remaining within the limits of the savanna biome. The latitudinal gradient follows a rainfall gradient with the driest savanna areas occurring in the lower latitudes with the northernmost site receiving ∼350 mm rainfall and the southernmost sites receiving ∼1,140 mm (**Table 1**).To account for the role of temperature and to extend the range of temperatures that the plants were exposed to, we paired each low altitude transplant site across the 850 km latitudinal gradient with a cooler higher elevation savanna site (above 500 m.a.s.l). Within each site, trees were planted with and without grass competition as grass competition is the strongest form of biotic competition in savannas (Riginos, 2009; February et al., 2013).

If climate is the primary limitation on the distribution of our tree species across a latitudinal gradient, then their distribution limits should coincide with climatic conditions where plant performance declines. In line with the literature, we predict that in the absence of grass competition, plant performance will decline with increasing distance from the range edge (Sexton et al., 2009; Hargreaves et al., 2013; Lee-Yaw et al., 2016). In the presence of grass competition, performance declines should match the current distribution limits and plants will not establish


at sites outside their range, a result that is frequently recorded in similar experiments in temperate systems (Hargreaves et al., 2013; Lee-Yaw et al., 2016). However if climate (either rainfall or temperature) is not the primary process limiting the current range we predict that the plants will grow and survive both inside and outside their distribution ranges in the absence of fire and herbivory. We excluded fire and herbivory at our experimental sites so as to test for climatic limits alone. However as the degree of fire and herbivory also change across rainfall and temperature gradients we used growth data from the experiment to develop models that link growth rates and observed fire return periods to explore how fire regime might interact with growth constraints to limit species range.

## MATERIALS AND METHODS

## Study Species

We transplanted two common savanna trees, Colophospermum mopane (mopane) (Kirk ex Benth.) Kirk ex J.Léon and Senegalia nigrescens (Oliv.). C. mopane is a highly dominant leguminous tree or multi-stemmed shrub, belonging to the detarioideae subfamily of the Fabaceae. Of the 1.5 million km<sup>2</sup> of the savannas in southern Africa, C. mopane is estimated to almost singly dominate 550,000 km<sup>2</sup> (25–35%) of savannas (Mapaure, 1994; Timberlake, 1995) where it often forms monospecific stands (Timberlake, 1995) on both nutrient rich and nutrient poor substrates. While dominating many thousands of hectares across southern Africa, its southern-most distribution ends abruptly in South Africa at ∼-22.5S in the west and at ∼24.0S in the east (Mapaure, 1994) (**Figure 1**). S. nigrescens (Fabacae) or the Knob thorn is a woody tree that can grow up to 20 m in height, although it usually remains between 8 and 10 m (**Figure 6**). It is most common and widespread in low-altitude savannas occurring from Tanzania to South Africa. The southern distribution limit for this species is at ∼28.5 degrees latitude and it is considered frost sensitive, preferring low-lying areas (Smith et al., 2012).

## Transplant Experiment

Eight transplant sites were established within the South African savanna biome across an 850 km latitudinal gradient (∼8 ◦ lat.), covering a rainfall gradient from 340 to 1,130 mm (**Table 1**). Four low elevation sites (<500 m.a.s.l) were paired with four high elevation sites (>500 m.a.s.l) at the same latitude. The paired cooler higher elevation sites were chosen to ensure adequate representation of cold and frost prone sites (Schulze, 2007). The latitudinal and altitudinal gradients included sites inside and outside the distribution ranges of both C. mopane and S. nigrescens (**Figure 1**). Sites across these gradients were selected if they occurred within a savanna, had a perennial, undisturbed grass layer and a clay content of 20–30% (both on the surface and at 20 cm). At each site, a 30 × 30 m plot was established and fenced to prevent mammal herbivory. Each plot was divided into 6 sub plots (**Figure 2**). Three replicated sub-plots were randomly cleared of all grass and forb biomass by the roots using a hoe whilst grass cover in the remaining three was preserved.

S. nigrescens and C. mopane seeds were collected from within their respective distribution ranges. S. nigrescens seeds were

collected within a ∼30 km radius of the South African Wildlife College (24.453 S, 31.405 E) and C. mopane seeds were collected within a 30 km radius of All days town (22.67 S, 29.10 E). The collected seeds were mixed and germinated at a nursery in the lowveld savanna at the South African Wildlife College. Seeds were planted in 1 l potting bags in a sand: clay mixture which was inoculated with soil taken from the roots of each species growing in their natural habitat. Following successful germination, the 480 one month old seedlings were transported to the sites and were planted in the early growing season (3 November−21 November 2010). Within each replicate five C. mopane and five S. nigrescens one month old seedlings were planted in random order one meter apart (**Figure 3**). To reduce losses to transplant stress and to allow the seedlings to establish, an equivalent of a 10 mm rainfall event was supplied every week for 1 month. Recruitment in these species is episodic (Botha, 2006; Stevens et al., 2014a) and by watering at the seedling establishment period

subplots of a cleared treatment, where all grass was removed. Five plants of each species were planted in each replicate, at one meter intervals.

we simulated an early growing season with good rainfall that was likely to result in successful establishment. As the requirements for seedling establishment have been documented we replicated these conditions so we could examine the determinants of plant success and growth following early establishment. For successful early seedling establishment the plants require warmth and moisture. These environmental conditions were met at all the savanna sites both inside and outside the distribution range of both species and were therefore not the limiting factor to establishment (Botha, 2006; Stevens et al., 2014a).

## Seedling Survival

We monitored plant growth for 2 years. Sites were revisited three times annually at the start, middle and end of the growing season (early November, January, April). A round trip between sites was a distance of 5,000 km; therefore more frequent measurements were not taken. During each measuring period seedling survival was assessed and the standing height and stem diameter of each seedling was measured. At the end of the growing season the grass biomass in the grass treatments was measured using a disc pasture meter (DPM). DPM readings were converted to biomass using the conversion developed by Trollope and Potgieter (1986). i-Buttons were placed at each site at 1.2 m above the ground and at 0.15 m above the ground to assess the thermal environment experienced by the seedlings. I-Button data was used to establish the site specific temperatures reported in **Table 1**. I-Buttons were suspended inside 4 mm thick PVC tubes (20 × 40 cm), and remained in the field for the duration of the experiment.

## Data Analysis

The data were analyzed in R. 3.5 (R Core Team, 2018). To examine the survival of seedlings during the course of the experiment we used the survfit function to create Kaplan-Meier survival curves for the full set of seedlings distinguishing between treatment (grass: no grass), distribution (inside: outside the species range) and species. Log-rank tests were used to compare

the groups. To examine if there was an effect of transplant shock on seedling survival we ran all the Kaplan-Meier survival curves both with the default alpha, and with an alpha of−10 which reduces the weighting of the first part of the dataset.

We used mixed-effects Cox proportional hazards models using the coxme (version 2.2.10) package (Therneau, 2018) in R to quantify the effects of treatment, distribution, mean annual temperature and mean annual precipitation and species as predictors of seedling survival, as well as the interaction between treatment and species. To account for the ties in the data (due to the limited number of sample dates) we used the "Efron" method for calculating ties. We ran the model for the full dataset, and for a dataset that excludes the transplant period (i.e., only data after day 130 were used in the model). We included both plot and site as random effects (plot nested within site). Candidate models were compared and ranked using Bayesian information criterion (BIC). If the best candidate models differed with a 1BIC <2 we used model averaging to determine the best model parameters.

In order to determine the relationship between final plant height and climate parameters, mean annual temperature (MAT) and precipitation (MAP) were calculated for each site using temperature data gathered from i-Buttons on site and rainfall records from either long term farm records or the closest SA weather service station. The final mean height differences between sites were analyzed using linear-mixed effects models in the nlme package (Pinheiro et al., 2009). Linear mixed effects modeling fit by maximum likelihood was used to assess the significance of mean annual temperature, mean annual rainfall, grass treatment (grass or no grass) and range (inside or outside of distribution range) on the height of saplings after 2 years of growth. We included mean annual temperature and mean annual precipitation as fixed covariates in the analysis and distribution (inside or outside range) and treatment (grass or no grass) as fixed categorical factors. Site was included as a random factor and plot was nested within site. All possible interactions were considered. The candidate models were compared and ranked using Bayesian information criterion (BIC). The second best candidate model for both species differed with a 1BIC >2 so model averaging was not used and the first best fit model with the lowest BIC was selected. We present the final parameter estimates, standard errors, and confidence intervals to demonstrate the effect size of the different parameters of the top two models.

## Modeling Escape From the Fire Trap

To model the probability of sapling escape from the fire trap we used mean maximum growth rates for the two fastest growing plants of each species from the grass treatment, for each site (Wakeling et al., 2011; Bond et al., 2012). Growth rates were calculated from the beginning of the second growing season to the beginning of the third. We did not include the first season in the growth rate calculations as growth in the first season is characteristic of seedlings which have different allocation patterns and hence different above ground growth rates than that of established saplings (Higgins et al., 2000). We used the growth model from Higgins et al. (2000) to model the height gain of savanna trees, and hence time to reach a fire proof height (∼3 m):

$$h = h\_{\mathcal{Y}^{-1}} + \left(1 - \frac{h\_{\mathcal{Y}^{-1}}}{h\_{\text{max}}}\right) \mathbf{g}\_{\text{s}}$$

where g<sup>s</sup> is the growth rate of stems (cm year−<sup>1</sup> ) (year 2), hmax is the maximum stem height (10 m), and hy−<sup>1</sup> is the stem height in the previous year. Starting stem heights were set at a starting height of 20 cm, the mean stem height after year one across all the transplant sites. We kept the starting height the same across sites so that responses could be restricted to the sapling stage, not the seedling establishment stage. The model was run for 50 years.

As the model produced an output of tree height per year, we could determine the time it took a tree at each site to reach the "escape height" of 3 m, given the maximum growth rate at a site (Wakeling et al., 2012). We used maximum growth rates as they are more likely to be ecologically relevant than mean growth rates as predictors of sapling release and successful recruitment into the adult stage (Wakeling et al., 2011; Bond et al., 2012). We plotted this value against the median fire return interval (FRI) for each experimental site. FRI was calculated based on fire return intervals for low-human density African savanna regions (Archibald, 2016). This was estimated by fitting a Weibull distribution (Johnson and Gutsell, 1994) to fire interval data derived from the MODIS burned area product (see Archibald et al., 2009; Archibald and Hempson, 2016 for methods) for ∼ 1,000 points from undisturbed areas.

## RESULTS

## Site Characteristics

Sites spanned across a latitudinal gradient (**Table 1**). With increasing latitude mean annual precipitation (MAP) increased from ∼400 to 1100 mm and mean annual temperature (MAT) decreased by ∼3 ◦C (**Table 1**).The differences in rainfall between high and low elevation sites were negligible except for the lower latitude northern sites (VN & PM) where the low elevation site received 120 mm more per year. High and low elevation sites differed in their mean annual temperatures with the low elevation sites having higher mean annual temperatures, warmer winters and mild winter minimums. VN, the low latitude higher elevation site was the exception, as it is had the highest MAT, however the mean monthly maximums were higher at the low elevation site. This site is also the driest site receiving a MAP of 344 mm. The wettest site in the low elevation, higher latitude (UK) received a MAP of 1,133 mm. The coldest site was the higher elevation KW site with a MAP of 18.9 and a July mean minimum of 0.68◦C (the coldest month). KW was also the only site that received regular frost events during winter.

## Seedling Survival

We used Kaplan-Meier survival curves to show how distribution, treatment and species impacted seedling survival during the experiment (**Figure 3**). The presence or absence of grass had the strongest impact on seedling survival with grasses reducing survival probability by ∼50% [X 2 (1, 457)= 31.0, p < 0.0001]. This effect still remained significant after accounting for transplant induced mortality at the beginning of the experiment [X 2 (1, 457)= 28.0, p < 0.0001]. Interestingly there was significantly lower survival inside the distribution range [X 2 (1, 457)= 34.0, p < 0.0001], but the bulk of this effect occurred during the initial transplant period (first 130 days). However distribution range was not significant when the early stages of the experiment were given less weight (alpha−10) [X 2 (1, 457)= 1.3, p =0.253]. There were no significant differences in the survival rates between species with both species showing similar responses to distribution range and grass treatment.

For the first dataset (all time steps included) mixed-effects Cox proportional hazard models with the most support (based on BIC and model inference) included grass treatment and distribution range. Although this model had a weight of 0.51 there were two other important models—one that included MAP, and one that included MAT (**Table 2A**). As there were three models where 1BIC <2, we averaged these three models. The model averaged co-efficients indicate that treatment, distribution range, MAT and MAP all have explanatory power. The terms with the highest weighting were treatment and distribution with a weighting of 1, whilst MAP and MAT each had a weighting of 0.23. The final averaged model indicates that distribution range (p < 0.01) and treatment (p < 0.001) are the only significant explanatory variables in the model. The log of the co-efficient values show that plots without grass had ∼2 times the survival rate as plots with grass, and plots inside the distribution range had about half the survival rates as plants outside the range (**Table 3A**).

To see whether the effect of distribution was limited to the early transplant period, we ran the same model but excluded all data before day 130. In this case the top 4 models all had very similar weights (0.10–0.11) with distribution range, MAP, MAT and treatment appearing as important terms (**Table 2B**). We therefore averaged the model results from models where 1BIC <2 (**Table 3B**). Treatment was the most important variable and received a weight of 1, followed by MAT (0.73), MAP (0.59), distribution range (0.42), and species (0.5). In the final model the grass treatment was the only significant variable (p < 0.0001) and individuals planted with grass were ∼three times as likely to die (**Table 3B**).

## Final Sapling Height

We examined how the final height was explained by distribution range, MAT, MAP, grass treatment and all possible interaction terms. We used BIC values and model inference to select the best fit model. The best fit final mixed linear model for both C. mopane and S. nigrescens included only the treatment variable. We did not use model averaging for the results as the next best fit model for both species had a BIC>2. Plants were significantly shorter at the end of the experiment if they were grown in the presence of grass (**Table 4**, **Figures 4**, **5**).

Contrary to predictions, frost did not cause sapling mortality. KW was the only site which experienced frost. Several frosts during the dry season in year 1 caused a decline in height for both species, (**Figures 4**, **5**). Frosted trees of both species were topkilled, but both species resprouted from the base in the following growing season. Resprouts were always multistemmed.

As grass is an important factor in sapling growth we investigated how it impacts sapling growth across the transplant experiment. Grass biomass had a significant effect on tree sapling growth within the sites and across the latitudinal gradient for both C. mopane (p < 0.01) and S. nigrescens (p < 0.001) (**Figure 6**). Sapling growth in the presence of grass for both species was highest where the grass biomass was lowest, and


TABLE 2 | Best models (where BIC <2) assessing the relationship between variables and seedling survival for (A) the entire experimental period and (B) the entire experimental period but with the first 130 days receiving a lower weighting (Alpha = −10).

*Models are ranked based on differences in the Bayesian Information Criteria (BIC). BIC weight is the weight of each model. Range, inside or outside distribution range; MAP, mean annual precipitation; MAT, mean annual temperature; Species, species of seedling (S. nigrescens or C. mopane); treat, presence or absence of grass.*

TABLE 3 | Final parameter estimates (β), standard errors (SE) and confidence intervals of model averaging based on top models for (A) analysis of entire experimental period and (B) analysis excluding first 130 days to account for effects of tranplant shock.


*Significant parameters are shown in bold. Range, inside or outside distribution range; MAP, mean annual precipitation; MAT, mean annual temperature; Species, species of seedling (S. nigrescens or C. mopane); treat, presence or absence of grass.*

sapling height progressively declined as grass biomass increased across the sites for both species (**Figure 6**).

## Modeled Tree Growth and Escape From the Fire Trap

Considering that seedling survival and performance were not well explained by distribution range we then looked at processes preventing sapling recruitment into adult size classes. Grass biomass not only reduces sapling performance (**Figure 6**), but also provides the fuel that determines fire return intervals (**Figure 7**). High grass productivity, linked to higher rainfall, permits frequent fires. We therefore predicted that a combination of slower growth rates and more frequent fire return intervals might be a significant factor preventing sapling recruitment. Species growing inside their distribution ranges should be more likely to reach fire escape height within the fire free interval common to their area. Our model indicates that C. mopane fails to escape from the fire trap when maximum growth rates are low (caused by high grass biomass) (**Figures 7**, **8**), or when fire return interval is short (**Table 1**, **Figure 7**). Model results (**Figure 8**) indicate that within C. mopane's distribution range

#### TABLE 4a | Mixed linear model outputs showing the top two models selected using BIC values for *C. mopane*.


*Bold terms indicate significance. Treatment, the presence or absence of grass, was the only factor that significantly affected tree growth. Confidence intervals are presented in square brackets.*

TABLE 4b | Mixed linear model outputs showing the top two models, selected using BIC values for *S. nigrescens*.


*Bold terms indicate significance. Treatment, the presence or absence of grass, was the only factor that significantly affected tree growth. Confidence intervals are presented in square brackets.*

sapling growth rates with grass are fast enough to result in trees escaping the fire trap. Outside the distribution range, although actual growth rates can be higher, fires are more frequent, FRI intervals decrease, and C. mopane is unlikely to reach escape height before the next fire. Likewise for S. nigrescens, the probability of successful escape from the fire trap only occurred within the existing distribution range of this species. At these sites maximum tree growth rate in the presence of grass can out-pace the fire return interval (**Figure 8**).

## DISCUSSION

The performance of two common savanna trees did not decline with increasing distance from the range edge. Instead plants still

survived and grew well hundreds of kilometers outside their distribution ranges. This is contrary to the generally observed pattern that species distribution limits are often niche limits (Lee-Yaw et al., 2016). In our study species, and perhaps more generally in common savanna plants, climate axes of the niche were poor predictors of species performance in and outside range limits. Plant performance instead was most strongly determined by grass biomass. Whilst saplings still survived and grew within the presence of grass, the poorer performance with high grass productivity reduced the probability of saplings escaping the fire trap between successive fires. Our results and modeling suggest an alternative framework for the determinants of species ranges.

Prediction of ranges for these and potentially other savanna species requires an understanding of the interactions between tree and grass productivity, climate and disturbance (top-kill) which determine the ability of saplings to escape the fire trap at any given site.

As the plants survived and grew well past their southern limits, a classic interpretation would be that both these species were dispersal-limited since suitable, but unoccupied, areas existed beyond their range boundaries (Hargreaves et al., 2013; Lee-Yaw et al., 2016). While dispersal limitation does occur, it is not common (Sexton et al., 2009; Hargreaves et al., 2013), and it is highly unlikely that these two widely distributed, dominant savanna species are dispersal limited. The fruits of C. mopane are well dispersed by wind and water (Styles and Skinner, 1997; Mlambo and Nyathi, 2004), and allozyme frequencies are similar across Southern African populations, indicating that an effective mechanism of gene flow occurs (Villoen et al., 2003). S. nigrescens retains closed pods on the tree until ripe, when gusting winds, associated with the onset of rains, rip the pods off the tree where they open, scattering the seeds widely (Miller, 1994; Grant and Thomas, 2001).

It is often documented that herbaceous competition strongly suppresses sapling and tree growth (Riginos, 2009; February et al., 2013) and likewise the presence of grasses in our study has a strong suppressive effect on seedling survival and sapling growth. It is through this interaction where climate probably has the strongest indirect influence on tree establishment as grass biomass shows a strong positive increase with increasing rainfall (O'connor et al., 2001). In the juvenile stages, tree

*S. nigrescens* and (B) *C. mopane*. Closed circles show a species was planted within its distribution range, open circles show the plant was grown outside its distribution range. The dashed diagonal line indicates the minimum fire return interval required for each time-to-escape-height value. If the points fall above the 1:1 line, into the shaded area, the plants have a high probability of escaping from the fire trap. Both species, if planted within their distribution range, have a reasonable probability of escaping the fire trap. Estimated growth rates are too slow for saplings to escape the prevailing fire regimes outside their distribution range.

competition with grasses is likely to be highest as tree and grass roots occupy the same soil layers (Cramer et al., 2007). We therefore hypothesize that the distribution range of both of these tree species is a product of the interaction between sapling growth rates, grass competition, and the fire return interval which depends on grass biomass for fuel. Here, climate plays an indirect role through its interactions with fire return interval and grass biomass (O'connor et al., 2001; Archibald et al., 2009). Grass biomass can set the maximum annual tree growth rate with growth rates decreasing with increases in grass biomass. The probability that any sapling will successfully escape the fire trap is dependent on its maximum growth rates (Wakeling et al., 2011; Bond et al., 2012). The sapling escape time (the time it takes for savanna tree to reach 3 m) needs to occur within fire free intervals. If a fire burns a sapling it will be top-killed and will have to start growth from a resprout. Our model therefore matches the time taken to reach escape height against the prevailing fire frequency of an area. We suggest that the probability of successful establishment is very low when the time to reach escape height is greater than the prevailing fire return interval. The chances of establishment increase as the times become more closely matched. The highest probability of establishment therefore occurs when escape height can be easily reached within the window of the prevailing fire free interval (Higgins et al., 2000).We therefore suggest that disturbance (e.g., fire) is the critical missing component in determining whether a plant can become a successfully recruiting adult, and hence should be explicitly considered when predicting the future ranges of species. Fire regimes can, and are, being modified by land users (Archibald et al., 2013). Whether climate change, or fire (and herbivore) management has greatest effects on changes in species distribution needs careful consideration in savannas.

Whilst we demonstrate our idea using fire return interval this framework applies to any process that causes top-kill in savannas including mammal herbivory (Staver et al., 2012) and frost (Holdo, 2007; Whitecross et al., 2012). In our experiment, frost caused top-kill in the regularly frosted site (KW). It is a climatic factor considered important in limiting the distribution of C. mopane (Henning and White, 1974; O'Connor and Bredenkamp, 1997; Burke, 2006). However we found that frost did not directly kill these seedlings, as often predicted, instead it acted as a powerful top-kill agent with a similar impact to that of fire. Both species experienced severe frost (**Figure 6**) and experienced complete topkill. In the subsequent growing season multiple shoots resprouted from the base and the plants continued to survive and grow but lost the previous season's height gains. In sites with regular frost, the equivalent metric would be the frost return interval, and saplings would need to reach a height threshold where frost damage is negligible (Whitecross et al., 2012) or rely on frost free intervals to ensure successful establishment (Holdo, 2007; Wakeling et al., 2012; Whitecross et al., 2012), or experience eventual population exclusion within that area (Wakeling et al., 2012). Likewise mammalian grazers can interact with the distribution of species through the removal of grass and the reduction of fire frequencies which in turn will lessen the fire imposed bottleneck (Hempson et al., 2017). Alternatively a high prevalence of browsers in the landscape has the potential to increase sapling top-kill further reducing the likelihood of escape from the fire trap (Staver et al., 2009).

Similar insights have been made at biome levels across continental scales where fire, vegetation and climate are recognized to be fundamental in predicting the distribution of the savanna biomes (Whittaker, 1970; Bond et al., 2005; Lehmann et al., 2014). As these processes interact to shape the distribution of a biome it should not be surprising that the dominant species within these biomes are shaped by similar forces. Several unrelated studies have demonstrated that fire is an important

## REFERENCES


environmental filter and determinant of species composition (Staver et al., 2012; Pausas and Ribeiro, 2017). For example, fire exclusion can drive a switch in species composition or a shift in the dominant functional types e.g., evergreen species (Trapnell, 1959; Plas et al., 2013; Charles-Dominique et al., 2017). However, we do not know of other studies showing how fire regime can restrict the regional distribution of savanna tree species. Our results indicate that top-down processes are important in explaining the distribution of these two savanna tree species and processes that cause top-kill like fire, frost and herbivory can act as strong environmental filters and thus be as important in shaping plant distributions as temperature in temperate systems. Climate-based species distribution models are widely used as a tool in predicting the future of African plant species and as a tool to shape climate change policy across the developing economies of Africa (Heubes et al., 2011; Bocksberger et al., 2016). We suggest that their assumptions need much wider testing given the growing evidence for consumer-control of these ecosystems.

## AUTHOR CONTRIBUTIONS

NS and WB conceptualized the idea. NS and SA designed the experiment and NS performed field work. SA commented on drafts, provided funding, fire data and advised on statistics, and performed some statistical analysis. WB provided funding, and commented on drafts.

## FUNDING

Funding was received from the Andrew Mellon foundation, NRF Thuthuka, CSIR Land Atmosphere feedback Parliamentary Grant and the fund for Spatial Planning for Protected Areas in Response to Climate Change (SPARC, Conservation International).

## ACKNOWLEDGMENTS

Alweet Hlungani, Henry Malan, Brian Blevin, Trevor Keith, Gail Bowers Winters and Prof Kevin Kirkman are thanked for their assistance and access to land. SANParks, Dinokeng and De Beers are thanked for access to their land. Thanks to Heath Beckett, Julia Wakeling and Chris Barichievy for assistance with the experiment. This submission has been adapted from a chapter in NS PhD thesis (Stevens, 2014).


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Stevens, Archibald and Bond. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Intraspecific Trait Variability in Andropogon gerardii, a Dominant Grass Species in the US Great Plains

Seton Bachle<sup>1</sup> \*, Daniel M. Griffith<sup>2</sup> and Jesse B. Nippert <sup>1</sup>

<sup>1</sup> Division of Biology, Kansas State University, Manhattan, KS, United States, <sup>2</sup> Forest Ecosystems and Society, Oregon State University, Corvallis, OR, United States

The climatic conditions in the North American Great Plains are highly variable, characteristic of an inter-continental climate. Antecedent climate history has impacted the flora of Great Plains grasslands, resulting in high species richness as well as dominance by only a few grass species, such as Andropogon gerardii. While the productivity of A. gerardii is well described, the individual physiological, and morphological characteristics that confer species dominance over wide spatial gradients are not clearly understood. We performed a literature search to assess intra-specific trait variability of A. gerardii from as many locations as possible. Ultimately, only 13 locations in the Great Plains have reported common plant functional traits (PFTs) for this species. To best represent site-specific climate conditions, plant functional trait data (8 PFTs) were collected from literature reporting ambient growing conditions, and excluded experimental manipulations. For most PFTs, we found insufficient data to fully quantify the range of variation across the geographical extent of A. gerardii dominance. This is surprising given that we focused on the most abundant grass in one of the most well-studied regions globally. Furthermore, trait data collected from our literature search showed a high degree of variability, but no strong relationships were observed between mean trait values and climate predictors. Our review of the literature on A. gerardii suggests a role for trait variability as a mechanism enabling the dominance of this species across large regions such as the Great Plains of North America.

### Edited by:

Oana Moldovan, Emil Racovita Institute of Speleology, Romania

### Reviewed by:

Cristian Atala, Pontificia Universidad Católica de Valparaíso, Chile Christine McAllister, Principia College, United States

> \*Correspondence: Seton Bachle sbachle@ksu.edu

#### Specialty section:

This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution

Received: 30 April 2018 Accepted: 03 December 2018 Published: 14 December 2018

#### Citation:

Bachle S, Griffith DM and Nippert JB (2018) Intraspecific Trait Variability in Andropogon gerardii, a Dominant Grass Species in the US Great Plains. Front. Ecol. Evol. 6:217. doi: 10.3389/fevo.2018.00217 Keywords: plant functional traits, grasslands, variability, ecophysiology, intraspecific trait variability

## INTRODUCTION

Biomes are often spatially delineated with sharp boundaries and attributed functional properties based on the primary vegetation represented (Bailey, 1998, 2005). These vegetation types are comprised of species that likely exhibit variation in functional traits that may respond to climatic gradients or change through time. Interspecific variation has been used as a proxy for understanding climate change, because it represents how plant functional types and ecosystem function may be altered (Ryan, 1991; Adler et al., 2006; Taylor et al., 2014). However, the focus on mean-trait values in the literature suggests that the true magnitude of intraspecific trait variation is commonly overlooked, and is rarely incorporated into climate and vegetation process models (Lambert et al., 2011; Johnson et al., 2015; Funk et al., 2017), although intraspecific trait variability appears promising for predicting species change (Lu et al., 2017). This is partly due to a lack of empirical work focusing on patterns of ecotypic and intraspecific trait variation (Violle et al., 2012). Intraspecific trait variation is the occurrence of a genotype expressing various phenotypes in a given environment; it's a combination of genetic (i.e., evolution) and environmental factors reflected in the variation within populations (forming ecotypes) (Valladares et al., 2007, 2014; Turcotte and Levine, 2016; Barbour et al., 2018). Measuring intraspecific trait variation allows for an indepth understanding of a species' ability to respond and adapt to environmental changes (Molina-Montenegro et al., 2018). In this review, we focus on the average trait expression (reflected in functional traits) change in a single species spanning a large continental gradient. Our goal is to assess the degree to which intraspecific trait variation in dominant species might contribute to the functional responses of grassland ecosystems.

The Central and Great Plains region of the United States experiences a continental climate of temperature extremes and both intra- and inter-annual variability in precipitation (Borchert, 1950; Weaver, 1968). A noteworthy characteristic of the grasslands in this region is the high floristic richness (Collins and Calabrese, 2012), yet dominance by a few C<sup>4</sup> grass species that encompass the majority of annual production (Dietrich and Smith, 2016). A combination of site-level and regional landscape heterogeneity contributes to genotypic and phenotypic diversity within dominant species (Olsen et al., 2013). This allows for a mosaic of different genotypes of dominant grass species to exist across the Great Plains region.

Previous research has shown the occurrence of broad genotypic and anatomical differences within species across regional gradients in the United States (Avolio and Smith, 2013; Olsen et al., 2013; McAllister et al., 2015). These differences may arise from climatic events causing changes in populations of dominant species (Hoover et al., 2014a, 2015; Hoffman et al., 2018). Functional trait variability may also play a large role in muting the negative impacts of large stress events, such as drought, across a landscape. Therefore, an innate advantage of increased resistance to stress events exists, and high functional trait variability may decrease the likelihood of species loss on a regional scale (Smith and Knapp, 1999).

The aim of this review is to measure within-region trait variation from a dominant species found in the Great Plains (temperate grassland; **Figure 1**). Focusing on several functional traits commonly used in the literature (specific leaf area, water potential, photosynthesis, stomatal conductance, aboveground productivity (individual and g/m<sup>2</sup> ), and below-ground biomass (individual and g/m<sup>2</sup> ), we have summarized the variability of functional traits within the species Andropogon gerardii Vitman. A. gerardii is a perennial C<sup>4</sup> grass that has significant genotypic and functional trait variability, facilitating a broad distribution throughout the Central and Great Plains. PFTs collected were analyzed to understand if the variability in a given trait varies/relates to gradients in multiple climate factors across this region, including mean annual precipitation (MAP), mean annual temperature (MAT), mean annual minimum temperature (MMinT), and mean annual maximum temperature (MMaxT) from 1980–2015. For example, we would expect biomass to vary across a climate gradient because productivity positively increases in this region with corresponding increases in precipitation and temperature (Nippert et al., 2006; Hufkens et al., 2016). Identifying trait variability within a single species may provide insight for the potential role of adaptive trait variability as a driver of population persistence across broad climatic space.

## LITERATURE SEARCH CRITERIA FOR Andropogon gerardii PFTS

We conducted a literature search to collect specific PFT data from a widely distributed grassland species: Andropogon gerardii, attempting to lend insight into one of many potential reasons for why some species achieve dominance over broad climatic space. A. gerardii was chosen because it is a dominant species in the tallgrass prairie and encompasses roughly 70% of the total aboveground biomass in grasslands throughout the Great and Central Plains region (Rogler, 1944; Weaver, 1968; Smith et al., 2017). We selected 6 important functional traits that reflect major axes of leaf economic variation and properties relevant to ecosystem function. The PFTs included in the literature search were specific leaf area (SLA), water potential (WPall; predawn and midday), photosynthesis (A), stomatal conductance (gs), above-ground biomass (AGB), and below-ground biomass (BGB). ABG and BGB data included biomass on an individual and a per square meter basis. Starting with SLA, we searched Google Scholar and Web of Science (which yielded identical results) with the string ["Andropogon gerardii" AND ("LMA" OR "SLA") AND "Great Plains"], which produced 60 search results. Only 12 of these studies (from 6 study sites) reported SLA. Similar searches were performed for the less common functional traits. SLA data was also collected from the TRY global database of plant traits (Kattge et al., 2011), which mirrored previous search parameters with greater success. This process resulted in a PFT dataset collected from 36 separate studies, including data from 17 research locations in 7 states from 1984–2017 (**Figure 1**). Our goal was to include data from across the U.S. Midwestern region and to exclude redundant data that were from the same projects, and only include "control" or "ambient" conditions. Data collection from literature varied by state. Kansas and Illinois contributed 59 data points (each point representing a single datum), which is a substantial amount of the total dataset. In order to determine differences in geographic location, statistical analysis was conducted via ANOVA before and after the data were normalized by natural log. MAP, MAT, MMinT, MMaxT (mean from 1980–2015), and geographic location of each data point was collected from PRISM (PRISM Climate Group, 2014) and used as fixed variables while each PFT was used as response variables. AICc model selection was also used to determine the most impactful climate parameters in the model using the "MuMIn" package (Barton, 2018), according to Grueber et al. (2011). Analyses were conducted in the statistical program R V3.4.3 (R Core Team, 2017). The geographic map (**Figure 1**) was produced using the "raster" package in program R (Hijmans, 2017).

## PLANT FUNCTIONAL TRAITS REFLECT ECOPHYSIOLOGICAL PROCESSES

Plant functional traits (PFTs) are commonly used to identify species' differences in growth, allocation, and competition in relation to environmental effects to reflect plant economics (Grime, 1979; Edwards et al., 2007; Guo et al., 2017; Volaire, 2018). PFTs represent morphological and physiological adaptations that often predict plant responses to biotic (competition, herbivory, etc.), and abiotic factors (MAP, MAT, etc.). PFT's typically include whole-plant traits, tissue specific traits (leaf, stem, and root), and physiological measurements (photosynthesis, stomatal conductance, transpiration, and water potential) (Pérez-Harguindeguy et al., 2013; Carmona et al., 2016). For instance, PFTs have been used to predict how individual A. gerardii and populations of other grassland species will respond to projected drought conditions (Chapin et al., 2000; Nippert et al., 2009; Volder et al., 2010; Liancourt et al., 2015; Maréchaux et al., 2015; Skelton et al., 2015; De La Riva et al., 2016).

The biological link between ecophysiology and environmental factors aids in predicting how species will respond to climatic changes (Nippert et al., 2011; Ocheltree et al., 2012; Hoover et al., 2014b; Griffith et al., 2016); therefore, such traits reflect the ability of a species, like A. gerardii, to respond to changing climatic conditions that are found in the Great Plains (Grime, 2001; McGill et al., 2006; Butterfield and Callaway, 2013; Losapio and Schob, 2017). Andropogon gerardii in the Great Plains exhibits many drought tolerant traits that allow either resistance and/or resilience in response to the abiotic stressors such as increased temperature and precipitation variability (Hoover et al., 2014b; Hoover and Rogers, 2016). Traits that are commonly correlated with drought tolerance include increased water-use efficiency (WUE), decreased leaf area (LA), higher specific root length (SRL), and lower turgor loss point (Eissenstat et al., 2000; Ripley et al., 2007; Hameed et al., 2012; Bartlett et al., 2014). Drought tolerant traits are likely not static in species with populations that span regional gradients, as there is adaptive benefit for greater trait variability in an environment that experiences high climate variability (Chapin, 1980; Avolio and Smith, 2013; Funk et al., 2017). For instance, WUE variability enables populations to maintain relatively high fitness with varying levels of water availability and mean temperatures (Briggs and Knapp, 2001; Nippert et al., 2007). Nippert et al. (2007) showed physiological trait variability for C<sup>4</sup> grasses provided an advantage for fast growth under favorable conditions and the ability to withstand (resist) drought during poor (stressful) conditions (Briggs and Knapp, 2001; Nippert et al., 2007). Trait variability may serve as a climate buffering mechanism (Valladares et al., 2007), which may be observable on the individual physiological scale (i.e., WUE), but also in other PFTs at the regional scale.

## TRAIT VARIATION IN Andropogon gerardii

Plant functional traits (PFTs) are known to differentiate between species, due to their evolutionary history (Violle et al., 2012; Cornwell et al., 2014; Valverde-Barrantes et al., 2017). However, less work has focused on intraspecific trait variability across large geographic scales. To emphasize this point, the dominant Great Plains grass species (Andropogon gerardii) is arguably the most well-studied dominant grass species in the Great Plains and specific leaf area (SLA) is the most widely reported functional trait, yet we only found 12 studies from 6 research locations where PFT data were reported from natural populations without experimental manipulation. In fact, three traits accounted for nearly half the data points in the study (**Figure 2A**), supporting the clear need for greater reporting of PFT data within common species that span large climate gradients (**Figure 2B**). Because single PFT data are under-represented in the literature, we aimed to incorporate multiple PFTs for A. gerardii to better understand intraspecific trait variability. More specifically, this amalgamation of PFTs may provide insight into the role of trait variability as a driver of plant species functioning over large regional scales. We assume that for a single species to extend over large geographic regions such as the Great Plains, the species would inevitably maintain a highly plastic phenotype at the population level which may buffer the whole species from variable climate conditions (**Figure 2B**). **Figure 2C** theorizes that populations from species with increased trait variability permits survival across a greater range of environmental conditions. For instance, if A. gerardii populations were subjected to drought conditions, varying levels of drought tolerance would be observed due to varying leaf water potential at turgor loss point (Maréchaux et al., 2015). Phenotypic variability has been observed to assuage effects from harsh abiotic pressures, however this may only be realized in the short-term (Becklin et al., 2016).

Long-term persistence of abiotic pressures will ultimately cause population reductions, due to climatic conditions moving out of the of historical climate parameters; this was observed in the droughts of the1930's (Romm, 2011; Becklin et al., 2016). Evidence from rain manipulation experiments provides insight into potential responses to climate extremes (Fay et al., 2002; Knapp et al., 2002; Nippert et al., 2009). For instance, Hoover et al. (2015), indicated that C<sup>4</sup> grasses subjected to drought conditions in the Colorado Plateau (35% reduction of annual rainfall) were observed to maintain cover for the first year, but decreased cover and increased mortality with prolonged exposure. Climate buffering can also be observed in a similar experimental design (rainfall manipulation) within the tallgrass prairie. A. gerardii did not display the same negative responses as other similar C<sup>4</sup> grasses to increasing climate variability, instead a relatively static response was observed (Fay et al., 2003; Avolio and Smith, 2013).

Increased trait variability within a species may provide physiological benefits in regions with high climate variability. A. gerardii lacks highly specific growing conditions and exists across broad geographic gradients in the U.S, including regions with hot and dry climate conditions. We hypothesize that A. gerardii and other generalists that dominate large geographic regions can be represented in **Figure 2B** as the blue line, whereas species that require more specific growing conditions are hypothetically represented as the black line. The theoretical curves in **Figure 1B** were created by using a low standard deviation (sd = 1; reflecting ∼68% of data explained) and a higher standard deviation (sd = 3; reflecting an increased data distribution) assuming an underlying normal distribution. Specialized species (black line) would not be capable of expanding over large heterogeneous landscapes due to the inability to withstand large fluctuations in temperature or precipitation (Linder et al., 2018), which is represented by a higher/narrower trait density (**Figure 2B**). Trait variability is documented to lead toward a more stable system due to niche stabilization which affects community composition, the function of the ecosystem, and response to abiotic factors (Turcotte and Levine, 2016).

The PFT data from Andropogon gerardii varied widely between research sites, but without discernable trends due to climate parameters like increasing precipitation and temperature (locations and PFT types are found in the **Supplemental Material**). For example, A. gerardii from Konza Prairie (KS), exhibited photosynthetic rates that included both the maximum and minimum of observed rates from all states included in the literature search, with many data points falling along the mean. Statistical results (ANOVA) show little to no discernable trend in PFTs when considering climate conditions at the geographic locations, which supports the concept of intraspecific trait variation allowing a single species to occupy such a large geographic range containing large precipitation and temperature differences. Photosynthetic rates did not vary by location, neither did stomatal conductance, SLA, or water potential (P > 0.05). Only two PFTs were observed to statistically vary by research location: above and below-ground biomass (P < 0.05). This result may reflect the small sample size from the literature as only 15 total data points were found for both individual and square meter collection methods. AICc model selection was performed on models containing all combinations of the PRISM climate variables to find the best model (given all variables) using the "model.sel" function, within the MuMIn package (delta <2, Royall's 1/8 rule, and cumulative sum of model weights were used to identify uncertainty and differences in the models). The full model including location, MAP, MAT, MMinT, and MMaxT best explained the SLA results gathered from the literature search (AICc = 974.6; weight = 0.939); which was 23 times more likely to be the best explanation for variation compared to the next model (AICc = 980.9; weight = 0.040) that did not include location as a parameter. Production above/belowground can exhibit a positive relationship (Nippert et al., 2006) with soil moisture and precipitation or relatively no change (Zhou et al., 2009), providing potential for variation across precipitation gradients within the Central and Great Plains. No trends due to climate parameters were visible in the analysis based on geographic location (**Figure 1**), but we speculate that additional production data collected across this regional gradient may result in a positive relationship with regional precipitation gradients similar to results from precipitation manipulation experiments (Fay et al., 2000, 2003).

The results from our literature review followed the ecological hypotheses presented in **Figure 2B**. PFTs were normalized (natural log) and combined to view large scale trends in the data. A. gerardii follows very closely to the "generalists" species parabola, meaning that PFTs are variable or plastic across a broad range of environmental conditions (**Figure 2C**). A. gerardii also expresses a low density of several PFTs, which indicates that traits are not static, but variable. Plastic PFTs allow A. gerardii the ability to respond positively in a given population, which buffers the A. gerardii species as a whole. This interpretation may contradict previously held claims that A. gerardii will experience geographic shifts or

experience large population reductions due to climate change (Gray et al., 2014; Smith et al., 2017). Thus, more research is required to understand the role that intraspecific variation plays in the expansion and survival of A. gerardii in the Great Plains.

## FUTURE DIRECTIONS

The literature containing plant functional traits covers many different types of ecosystems and hundreds of species (Pérez-Harguindeguy et al., 2013), which have been emphasized in climate change literature encompassing major biomes across the world (Liancourt et al., 2015). Here, we used PFTs to identify inherent trait variability within a single species, and identify a potential role of adaptive variability as a driver of species persistence across a regional climate gradient. Results from this literature review suggest that the PFTs observed in A. gerardii do not statistically differ (excluding AGB and BGB) between locations measured (**Figure 1**). Findings from this review underline the importance of adaptive trait variability to permit greater phenotypic plasticity, which provides population buffering for some species that exist across broad climate gradients.

Moving forward, an increased focus from interspecific to intraspecific species trait variation may provide a greater understanding of how future climate variability will impact native plant species that span large regional scales. More specifically, A. gerardii is the quintessential prairie species; yet an extensive examination of the literature showed a relatively small number of sites reporting trait data. We advocate for the development of a grass trait network to examine the effect of climate change on specific dominant species in grasslands worldwide. This network could follow the framework created by Nutrient Network (NutNet; http://www. nutnet.org) and Drought Network (DroughtNet; http://www. drought-net.colostate.edu) to standardize measurements and

procedures, and allow for more consistent interpretation of PFTs response to changes in abiotic factors. Increased documentation of spatial climate gradients and species distributions will increase our understanding of the role of trait variability in species resistance and resilience to future changes in the environment.

## AUTHOR CONTRIBUTIONS

SB and DG were both responsible for collecting data for the literature search. SB was responsible for writing the manuscript and collecting the bulk of the data. DG gave insight into the manuscript and utilized the data to create the precipitation map. JN served as SB's advisor,

## REFERENCES


aided in the literature search, and helped shape the manuscript.

## ACKNOWLEDGMENTS

We would like to thank the Konza Prairie LTER program (NSF DEB-1440484) and the NSF Dimensions of Biodiversity program (NSF 002893) for funding this research.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2018.00217/full#supplementary-material


along a large-scale aridity gradient in northeastern China. Sci. Rep. 7:40900. doi: 10.1038/srep40900


of rapid adaptive evolution in seed traits? Evidence from a latitudinal rainfall gradient. Front. Plant Sci. 9:208. doi: 10.3389/fpls.2018.00208


mycorrhizal association on the functional trait variation of fine-root tissues in seed plants. New Phytol. doi: 10.1111/nph.14571


Zhou, X., Talley, M., and Luo, Y. (2009). Biomass, litter, and soil respiration along a precipitation gradient in southern great plains, USA. Ecosystems 12, 1369–1380. doi: 10.1007/s10021-009-9296-7

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Bachle, Griffith and Nippert. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.