Woody plant encroachment into grasslands: spatial patterns of functional group distribution and community development.

Woody plant encroachment into grasslands has been globally widespread. The woody species invading grasslands represent a variety of contrasting plant functional groups and growth forms. Are some woody plant functional types (PFTs) better suited to invade grasslands than others? To what extent do local patterns of distribution and abundance of woody PFTs invading grasslands reflect intrinsic topoedaphic properties versus plant-induced changes in soil properties? We addressed these questions in the Southern Great Plains, United States at a subtropical grassland known to have been encroached upon by woody species over the past 50-100 years. A total of 20 woody species (9 tree-statured; 11 shrub-statured) were encountered along a transect extending from an upland into a playa basin. About half of the encroaching woody plants were potential N2-fixers (55% of species), but they contributed only 7% to 16 % of the total basal area. Most species and the PFTs they represent were ubiquitously distributed along the topoedaphic gradient, but with varying abundances. Overstory-understory comparisons suggest that while future species composition of these woody communities is likely to change, PFT composition is not. Canonical correspondence analysis (CCA) ordination and variance partitioning (Partial CCA) indicated that woody species and PFT composition in developing woody communities was primarily influenced by intrinsic landscape location variables (e.g., soil texture) and secondarily by plant-induced changes in soil organic carbon and total nitrogen content. The ubiquitous distribution of species and PFTs suggests that woody plants are generally well-suited to a broad range of grassland topoedaphic settings. However, here we only examined categorical and non-quantitative functional traits. Although intrinsic soil properties exerted more control over the floristics of grassland-to-woodland succession did plant modifications of soil carbon and nitrogen concentrations, the latter are likely to influence productivity and nutrient cycling and may, over longer time-frames, feed back to influence PFT distributions.

Introduction

Increases in the abundance of trees and shrubs have been reported in ecosystems world-wide [1-3]. The resultant change in physiognomy from grassland, steppe or savanna to shrubland or woodland has significant impacts on ecosystem productivity, trophic structure, nutrient cycling and biodiversity [4-8], and may have significant implications for local [9], regional [10,11] and global [12] carbon cycling given the geographic extent of these ecosystems on Earth. Woody plant proliferation in recent decades has been reported in arctic tundra, in temperate, subtropical, tropical, coastal and montane grasslands, in hot and cold desert grasslands, and in savannas and steppe [3,4,13,14]. Because woody plant encroachment is a worldwide phenomenon, it is important to understand it in general terms if we are to anticipate and predict where and how the abundance of shrubs and trees might change under current and future environmental conditions. To date, although there are studies focused on multiple invaded woody species [15-17], the majority of studies related to woody plant encroachment have focused on key encroaching species [11]. From such studies we can only speculate as to i) which of the many traits specific to those species make them so successful; and ii) the relative importance of those traits.

A qualitative survey of published literature indicates that encroaching woody species represent a broad range of plant functional types (PFTs), ranging from shrub to tree in stature; from evergreen to deciduous, malacophyllous to sclerophyllous, and broad-leaved to needle-leaved in leaf habit, structure and size; from N2-fixing to non-fixing; from deep rooted to shallow rooted; and from mesophytic to xerophytic in water relations [18-22]. However, robust generalizations regarding traits or suites of traits that may make some species and the PFTs they represent more successful than others in invading grasslands do not exist. Are some woody plant growth forms and PFTs better suited than others to invade grasslands; and if so, under what environmental conditions? Furthermore, it is well-known that woody plants may modify soils and microclimate subsequent to their establishment. To what extent does this induced change affect woody PFT distribution and abundance?

The Southern Great Plains of North America is a region where the conversion of grass to woody plant dominance over the past 100+ years has been well-documented by a variety of sources including diaries of early settlers, time-series remote sensing, tree rings, carbon isotopes, and ecosystem models [23]. Potential drivers of this shift in life form abundances include livestock grazing, fire suppression, climate change, and atmospheric CO2 enrichment [24]. Topoedaphic properties [25] and fluctuations in native herbivore abundance [26] may locally constrain or mitigate responses to these drivers. While the phenomenon of woody plant encroachment and its effects on soil C and N cycles has been widely examined in this region [9], details pertaining to the development of woody plant communities in contrasting topoedaphic settings are largely unknown. The woody flora of the subtropical Tamaulipan Biotic Province, which encompasses southern Texas and northern Mexico [27], is highly diverse [18] and includes an array of PFTs [28]. This high floristic diversity of woody vegetation affords the opportunity to determine if certain woody PFTs might be better suited than others for encroachment into grasslands. Here, we sought to determine similarities and differences in the PFT composition of four shrub communities with varying establishment time developing on former grasslands along a landscape-scale catena (hill-slope) topoedaphic gradient (Figure 1). Specifically, we asked: (1) What are the similarities and differences in species and woody PFT composition of these communities? (2) Do PFT overstory-understory relationships change along the gradient? and (3) What is the influence of intrinsic soil physical properties (e.g. texture) relative to that of plant-induced changes in soil chemical properties (e.g., organic carbon and nitrogen content) on PFT abundance and distribution? Because our study was conducted on a local scale we were able to control for climate and land use history, and factors that might otherwise confound comparisons of species and functional types.

Figure 1

Aerial view of the study site showing a transect spanning shrub cluster and grove communities in savanna parkland uplands, a woodland community in an intermittent drainage, and a playa savanna community.

Each of these communities is developing on former grassland. White dots are GPS locations at 5 m intervals. Dark and light gray colors indicate woody and herbaceous vegetation, respectively. Trees and shrubs were inventoried in contiguous plots (6 x 6 m and 2 x 2 m, respectively) along the transect.

Methods

Study site

Field studies were conducted at the Texas AgriLife La Copita Research Area (LCRA, 27°40’N; 98°12’W ), approximately 24 km southwest of Alice, Texas, USA. LCRA is owned by Texas A&M University and no specific permission is required to do research at this site. This study did not involve any endangered or protected species. Climate of the region is subtropical with a mean annual temperature of 22.4°C and mean annual precipitation of 680 mm. Annual precipitation is bi-modally distributed with maxima in May/June and September. Elevation ranges from 75 to 90 m. The site was a working cattle ranch with a long history of livestock grazing before its designation as a research area in early 1980s. Fire has been suppressed while it was a working ranch and later as a research site. In the 1980s and 1990s, various brush management practices (mechanical and herbicidal) were implemented on the site. The woody plant flora at this site consists of 49 species [19] encompassing a range of stature and ‘woodiness’ that ranges from trees (arboreal) to shrubs (fruticose) to subshrubs (suffruticose/suffrutescent); and a wide range of leaf traits (texture, longevity, specific leaf area, N content) and rooting depths [20,21]. Gas exchange, water relations and isotopic composition of the various plant functional types (PFTs) have been well-described [29-33].

Uplands with sandy loam surface soils (Typic Argiustolls) are characterized by savanna parklands (Fig. 1). A well-developed argillic horizon ~ 40 cm below the surface is laterally extensive and supports a grassy matrix containing scattered, discrete shrub clusters. Cambic inclusions with poorly expressed argillic horizons are dispersed throughout the uplands and support large groves of woody vegetation [1]. Convex uplands grade gently (1-3% slopes) into lower-lying intermittent drainages with clay loams and clays (Pachic Argiustolls) characterized by closed-canopy woodlands. Oval-shaped playas (Ustic Epiaquerts and Vertic Argiaquolls) with no external drainage occupy the lowest portions of the landscape, with vegetation ranging from grassland to open woodland [34,35]. Prosopis glandulosa var. glandulosa (hereafter “Prosopis”) is the dominant overstory in all woody communities. Zanthoxylum fagara , Condalia hookerii, Diospyros texana, Celtis pallida, and Berberis trifoliolata are common shrub species (plant nomenclature follows [36]). Upland grasslands consist primarily of C4 grasses in the genera of Aristida, Bouteloua, Cenchrus, Chloris, and Setaria; and a diversity of herbaceous dicots. Playa communities are typically characterized by a dense and continuous ground layer of the C4 grasses Paspalum pubiflorum var. pubiflorum and Bothriochloa ischaemum [34].

Field sampling and lab analyses

A 309-m transect traversing the major plant communities (grassland, discrete cluster, grove, drainage woodland, and playa) was established along the topoedaphic gradient extending from the convex upland downslope into the concave playa basin (Figure 1). Elevation was determined based on a topographic field survey conducted in April 2004.

Soils cores (15 cm deep, 2.24 cm diameter) were collected at 1 m intervals (2 cores per point) along the transect. One core was used to quantify bulk density (BD) and soil texture, and the other was used to determine soil organic carbon (SOC), total nitrogen (TN), and pH. Another set of soil samples was collected at 3-m intervals during a 2-day period (no precipitation either day) in April 2005 to determine soil volumetric water content (VWC). We recognize that the soil VWC is highly dynamic; and that a one-time measurement of soil VWC at a single time and shallow depth is highly superficial. Thus, we used VWC only as a rough indicator of potential spatial variation in moisture along the catena gradient. Woody species with few and relatively large (basal diameter > 5 cm) stems and a height potential of > 4 m were recorded as trees. Their basal diameters were measured in contiguous 6 m × 6 m cells (plots) centered along the transect (n=51). The size and abundance of shrubs with relatively small stems (basal diameter < 5 cm) was recorded within 2 m × 2 m cells (plots) centered on the soil sample points (n = 154) along the transect. The shrub category thus included juveniles of woody species with potentially arborescent stature. Species characterized by suffrutescent growth habits (herbaceous perennials with woody bases) were not included.

We characterized woody species into different functional types according to available references and field observations. N2-fixing PFT data is from Zizter et al. 1996 [22]; Leaf texture data is from Nelson et al. 2002 [21]; root depth information is from Watts 1993 [20] and Boutton et al. 1999 [37].

Soil VWC was determined by weighing soil cores before and after oven drying at 105°C for 24 hours. Bulk density and soil texture were quantified using the core method and pipet method respectively [38]. Soil cores designated for SOC and TN determination were dried at 60°C for at least 48 hours, passed through a 2 mm screen to remove coarse organic fragments and gravel, and then pulverized to a fine powder in a centrifugal mill (Angstrom, Inc., Belleville, MI, USA). Samples were weighed into silver capsules (5 x 7 mm) using a microbalance, treated with HCl vapor in a desiccator to volatilize carbonate-carbon, dried thoroughly, and sealed into the capsules. SOC and TN concentrations were determined by dry combustion using an elemental analyzer (Carlo Erba EA-1108, CE Elantech, Lakewood, NJ). The concentrations were then converted to densities (g m-2) to a depth of 15 cm by multiplying by BD. Soil pH (Accumet Basic pH meter, Fisher Scientific) was determined using a CaCl2 solution (0.01 M CaCl2) containing 12 g soil.

Statistical analyses

Soil variables in each plant community type were compared using ANOVA with Tukey’s corrections for post hoc comparisons. Shrub and tree distributions along the topoedaphic gradient were explored using canonical correspondence analysis (CCA; [39]) in PC-ORD (version 5, MjM software, Gleneden Beach, OR; [40]). CCA was performed on two matrices: the plot by species matrix as the dependent variables, and the plot by environmental variables matrix as the independent (explanatory) variables. Environmental variables included soil texture (% sand, % clay), SOC, TN, pH, and BD. Silt was excluded because it is a linear combination of sand and clay (silt=100 – sand – clay), and VWC was excluded owing to our low sampling frequency. Matrices of total basal area per plot were analyzed separately for trees and shrubs because their plot sizes differed (6 m × 6 m and 2 m × 2 m, respectively). Environmental variables for each plot in the matrix were averages of their respective within-plot replicate measurements.

Partial CCA in PC-ORD was used to partition the explained variation of woody species composition and to quantify the relative influence of environmental variable groups on community structure. Inertia of ordination, which quantifies variation in species composition, is additive and can be distributed into groups of environmental variables [39]. Here, we used the proportion of total variation explained (TVE) by groups of environmental variables instead to quantify their relative contributions [41,42]. The relative influence of two groups of environmental variables was evaluated for the tree and the shrub species matrices. Prior to partial CCA, variables in each group were evaluated for their independent and significant contribution to variation in species composition [43] using randomization tests. Only significant (p < 0.05) variables were included in variation partitioning, as described in [44,45] (pH was excluded as it was non-significant). One group of environmental variables was comprised of soil texture (% sand, % clay) reflecting long-term pedogenesis along the hill-slope gradient. The second group of variables included SOC and TN, both of which are known to increase with time of site occupation by woody plants at this site [9,23]. We did not include BD in either group because it represents both long-term pedogenesis and shorter-term effects of organic matter additions and was not a strong explanatory variable

Results

Species composition and soil properties

Changes in elevation and edaphic properties along the hill-slope transect are depicted in Figure 2 and summarized in Table 1. On the date of sampling, soil VWC was highest in playa basins and the woodlands near these basins, and was uniformly lower elsewhere. Soil (0-15 cm) sand content was the greatest in the uplands, with clay content increasing down-slope along the catena and peaking in the playa basin. BD was highly variable along the hillslope. SOC density was also variable, but generally higher in intermittent drainage woodland and playa landscape locations and lower in uplands. Within uplands, the grassland community had uniformly low SOC with distinctive peaks in SOC being associated with grove and shrub cluster communities embedded within the grassland matrix. Spatial patterns of TN mimicked those of SOC (data not shown; correlation between SOC and TN =0.93; p<0.01).

Figure 2

Topoedaphic variations based on sampling at 1-m intervals along a transect extending from plant communities in a savanna parkland upland (grassland, shrub cluster, grove), through intermittent drainage woodlands and into a playa savanna community.

VWC = soil volumetric water content; SOC = soil organic carbon.

Table 1

Mean (+ standard error) values of soil (0-15 cm depth) variables and ANOVA results along a catena gradient in a subtropical grassland undergoing woody plant encroachment.

	Upland Grassland	Upland Cluster	Upland Grove	Drainage Woodland	Playa	ANOVAp
Bulk density (g cm-3)	1.4±0.01a	1.3±0.02ab	1.2±0.02b	1.2±0.01b	1.8±0.03b	<0.001
VWC (%)	5.1±0.1a	4.6±0.4a	4.2±0.1a	7.4±0.2b	16.4±0.5c	<0.001
PH	7.1±0.1a	7.1±0.2a	7.2±0.1a	7.2±0.1a	6.0±0.1b	<0.001
Organic carbon (g m-2)	1233±22a	1482±69ab	1925±87b	3035±99c	3142±155c	<0.001
Total nitrogen (g m-2)	123±2a	143±5ab	182±8b	288±10c	253±14c	<0.001
Sand (%)	79±0.1a	78±0.4a	78±0.4a	70±0.7b	48±0.8c	<0.001
Silt (%)	9±0.2a	10±0.3a	9±0.2a	13±0.2b	20±0.3c	<0.001
Clay (%)	12±0.1a	12±0.2a	13±0.2a	17±0.5b	32±0.6c	<0.001

Different superscripts indicate significant difference among communities. Number of samples: grassland=66, cluster=18, grove=59, woodland=125, and playa=41.

A total of 9 species with tree-stature potential occurred along the topoedaphic gradient (Table 2). The number of tree species encountered was lower in shrub cluster (4) and playa (6) communities; and higher in grove (7) and drainage woodland (8) communities. Prosopis trees in shrub cluster communities were dead, so were not included in these tallies. Among tree species, Acacia farnesiana (deciduous, potential N2-fixer; note: ‘potential’ is hereafter implicit in all uses of ‘N2-fixer’) and Z. fagara (evergreen, non-fixer) occurred in all four woody plant communities. Acacia rigidula (deciduous; N2-fixer) was restricted to woodland communities; and Parkinsonia aculeata (deciduous; no known report of N2-fixation, thus treated as non-fixer) was restricted to playa communities. Z. fagara had the highest stem density (183 ha-1), with four other tree species having stem densities >100 ha-1. Prosopis (deciduous, N2-fixer) had the greatest average individual plant stem basal area (543 cm2/plant).

Table 2

Occurrence of tree-stature species (> 5 cm basal diameter) within 6 × 6m contiguous plots along a transect spanning shrub cluster and grove communities in savanna parkland uplands, woodlands of intermittent drainages, and a playa (Figure 1).

Species	Cluster		Grove		Woodland		Playa
	N	BA	N	BA	N	BA	N	BA
Acacia farnesiana	2	0.6±0.1	2	0.2±1.4	8	1.6±3.9	19	17.0±3.3
Acacia rigidula	0	-	0	-	2	0.5±4.6	0	-
Celtis pallida	0	-	5	0.7±0.7	37	3.8±1.0	22	6.6±1.1
Condalia hookeri	6	2.5±1.3	5	0.7±2.0	39	2.4±0.6	0	-
Diospyros texana	0	-	3	0.3±1.1	50	2.1±0.3	5	0.5±0.2
Karwinskia humboltiana	2	0.6±1.8	4	0.3±0.3	1	0.1	0	-
Parkinsonia aculeata	0	-	0	-	0	-	1	0.46
Prosopis glandulosa	0	-	22	9.6±1.5	34	13.2±1.9	4	12.4±9.8
Zanthoxylum fagara	2	4.2±11.7	21	4.2±1.4	37	2.7±0.8	8	1.3±1.5

N = number of plants; BA = basal area (cm2 m-2; mean ± standard error). Prosopis glandulosa var. glanduloas was present in shrub clusters, but had died and hence was not recorded. See Table 3 for growth form and functional group categories represented by the species.

Table 3

Occurrence (N) and mean basal area (BA, mean + SE; cm2 m-2 ) of shrub-stature species within a 2-m wide belt transect (Figure 1).

Species	Functional attributes			Cluster		Grove		Woodland		Playa
	Leaves	N-fix	Roots	N	BA(cm2)	N	BA(cm2)	N	BA(cm2)	N	BA(cm2)
Acacia farnesiana$	D, M	*,#	Deep	0	-	0	-	2	30.8±1.4	1	22.1
Acacia greggii$	D, C	*,##	Int	5	30.7±0.3	0	-	0	-	0	-
Acacia rigidula$	D, C	*,#	Int	0	-	0	-	13	111.4±0.5	0	-
Aloysia gratissima$	D, M		Shallow	0	-	0	-	257	173.5±0.0	0	-
Berberis trifoliolata§¶	E, S		Shallow	19	31.3±0.2	10	10.9±0.2	20	17.5±0.1	0	-
Bernardia myricaefolia	D, M		-	0	-	0	-	1	0.5	0	-
Celtis pallida§¶	D, C	*,##	Int	6	1.8±0.1	109	530.1±0.1	562	1732.8±0.0	5	96.6±0.7
Colubrina texensis$	D, C	*,##	Deep	289	106.7±0.0	313	151.9±0.0	30	22.2±0.1	0	-
Condalia hookeri§$	FE,M	*,##	Int	15	73.4±0.5	18	9.7±0.0	12	187.4±0.7	0	-
Diospyros texana§$	FE,C		Int	0	-	12	28.0±0.3	136	854.6±0.1	7	115.5±1
Ephedra antisyphilitica	D,PS		-	0	-	0	-	1	2.8	0	-
Eysenhardtia texana	D,M	*,#	-	22	50.4±0.2	0	-	0	-	0	-
Forestiera angustifolia	E,M		-	2	3.5±0.3	1	2.8	4	9.1±0.2	0	-
Gymnosperma spp.	D,C		-	0	-	10	4.3±0.1	1	2.3	0	-
Karwinskia humboltiana	E,C	*,##	-	0	-	113	363.1±0.1	184	190.4±0	0	-
Parkinsonia aculeata	D, M,PS	*	Shallow	0	-	0	-	0	-	2	13.2±0.1
Prosopis glandulosa§¶	D,C	*,#	Deep	1	12.6	9	178.9±0.7	10	123.4±0.7	0	-
Schaefferia cuneifolia¶	E,S		Int	20	31.2±0.1	84	58.8±0	160	101.5±0	0	-
Zanthoxylum fagara§¶	E,C		Shallow	85	190.1±0.1	169	574.3±0.1	140	490.6±0.1	3	10.2±0.7
Ziziphus obtusifolia§¶	SD,M	*,##	Int	0	-	3	14.4±0.6	4	11.3±0.4	5	19.6±0.4

Functional attribute codes represent leaf habit, leaf texture, rooting depth and N2-fixation potential respectively. Codes are as follows: Leaf habit [E = evergreen, FE = facultative evergreen (leaves maintained during mild winters), D= deciduous, SD=summer deciduous, PS =photosynthetic stem]. Leaf texture [S = sclerophyllous (thick and stiff),C= coriaceous (thick and leathery), M = malacophyllous (thin and pliable)] (based on authors’ field observations or, where noted, Nelson et al. 2002)]. N-fixation [* = potentially N2-fixing (based on family affinities; e.g. Fabaceae, Rhamnaceae). and , = species for which N2-fixation potential has been assessed (Zitzer et al. 1996) and for which nodulation was found to occur (#) or not occur (##)]. Root depth [int = intermediate, Deep=deep, Shallow=shallow]. §,¶ ,$ = species for which information on leaf attributes (longevity, specific leaf area, N content, etc.; § = from Nelson et al. 2002) and relative rooting depth (¶ = from Watts1993; $ = from Boutton et al. 1999) are known.

A total of 20 shrub species were encountered within the 2-m wide belt centered on the transect line. All species with tree-stature potential (Table 2) also occurred as shrubs (Table 3). Individuals of A. rigidula (deciduous, N2-fixer), Aloysia gratissima (drought deciduous, non-fixer), Bernardia myricaefolia (deciduous, non-fixer) and Epheda antisyphilitica (stem photosynthesis, non-fixer) were encountered only in the woodland community; A. greggii (deciduous, N2-fixer) occurred only in the discrete cluster community; and shrub-stature Parkinsonia aculeata plants (deciduous, non-fixer) occurred only in the playa. In contrast, Z. fagara and C. pallida shrubs occurred in all four plant communities. Prosopis, D. texana, P. aculeata, and Acacia spp. were among the largest shrubs (average basal areas > 6 cm2 per stem), and C. texensis was among the smallest (average basal area per stem < 0.5 cm2).

Prosopis (deciduous, N2-fixer) was the dominant overstory plant in all woody plant communities (although it had died in the shrub clusters). C. texensis (deciduous, non-fixer) and Z. fagara (evergreen, non-fixer) dominated the shrub component in upland clusters and groves, whereas C. pallida (deciduous, non-fixer) and D. texana (deciduous, non-fixer) dominated the shrub understory in intermittent drainage and playa communities (Table 3).

Woody plant species richness (maximum encountered) was lowest in playa communities (n=7), highest in woodlands (n=17) and intermediate in shrub clusters and groves (n= 12 and 13, respectively) (Table 4A). Richness of deciduous species in lowland woodland and playa communities (n=9 and 5, respectively) exceeded that of evergreen species (n=5 and 1, respectively), whereas the richness of deciduous and evergreen species in upland cluster and grove communities was comparable (n=5 or 6). From a taxonomic perspective, encroaching woody plants encountered in our catena-scale sampling were primarily deciduous (60% of species); but evergreen (25%) and facultative evergreen species (10%) were also well-represented. About half of the encroaching woody plants were potential N-fixers (55% of species). However, only 4 species (20%) were confirmed to have nodulation in the field. Most potential N-fixing species were deciduous; and all deciduous species had either coriaceous (6 species) or malacophyllous (6 species) leaves. Leaf texture of evergreen species (n=7) ranged from coriaceous (n=3) to malacophyllous (n=2) to sclerophyllous (n=2). Rooting depths ranged from shallow (4 species) to intermediate (7 species) to deep (3 species); but were unknown for 6 species. Rooting depths for confirmed N-fixing species were either intermediate (1 species) or deep (2 species) (and unknown for the fourth species).

Table 4

Woody species functional groups and their relative basal area along the transect.

	A. Woody Species Richness						B. Woody species relative basal area (%)
Functional Traits	Cluster (n=12)	Grove (n=13)	Woodland (n=17)	Playa (n=7)	Total (n=20)		Cluster (540)	Grove (1943)	Woodland (4,089)	Playa (315)	Total (6,987)
Seasonal Leaf Habit [ a , b ]
Deciduous	6	6	9	5	13		38	46	55	60	51
Evergreen	5	5	5	1	5		48	52	20	4	31
Facultative Evergreen	1	2	2	1	2		14	2	26	37	19
Leaf Texture [ c ]
Malacohyllous	4	4	6	3[ c ]	7		24	1	10	23	10
Coriaceous	6	8	8	4	9		64	95	87	77	87
Sclerophyllous	2	1	2	0	2		12	4	3	0	4
Photosynthetic stem	0	0	1	1[ b ]	2		0	0	0	0	0
N-fixation [ d ]
N-Fixers	4	2	3	2	4		12	10	7	16	9
non N-fixers	8	11	14	5	16		88	90	93	84	91
Rooting Depth
Shallow	2	2	3	2	4		42	30	17	8	22
Intermediate	4	5	6	3	7		26	33	74	76	58
Deep	3	3	3	2	3		22	18	5	16	10
Unknown	3	3	5	0	6		10	19	5	0	9

(A) Number of woody species representing various functional groups in four tree-shrub communities developing along a catena gradient in a former grassland; and (B) their relative contribution (%) to total basal area in the community (values under column headings are total basal area (cm2 m-2) in the community. See Tables 3 for species and their functional trait attributes.

a Ephedra antisyphilitica (stem photosynthesis; woodland only) not included in count.

Summer deciduous species (one only; Zizyphus obtusifolia) included with ‘deciduous’

Parkinsonia aculeata has photosynthetic stems and malacophyllous leaves, so is counted in each category

Parkinsonia aculeata is not counted as a N-fixer

From a relative basal area (BA) perspective, woody communities developing on grassland were co-dominated by evergreen (including facultative evergreen) and deciduous PFTs (Table 4). Facultative evergreen PFTs were most abundant in the two upland communities. Species with coriaceous leaves dominated all woody communities along the catena gradient (> 64% of total BA) and were most strongly expressed in upland grove and lowland woodland communities (95 and 87% of total BA, respectively). N-fixing species were clearly in the minority, contributing only 7% (woodlands) to 16 % (playas) of the total BA. Basal area in upland clusters and groves was fairly evenly partitioned among shallow, intermediate and deep-rooted PFTs, whereas PFTs with intermediate depths dominated lowland woodland and playa communities (74 and 76% of BA, respectively).

CCA ordination

CCA ordination Axes I and II accounted for 18.0% and 4.4% of the total tree species inertia (3.42), respectively. The first CCA axis was strongly and positively correlated with soil clay, and negatively correlated with percent sand and pH. Tree species (Figure 3A) and communities (Figure 3B) were well-separated along this axis. There was considerable overlap among upland cluster, upland grove and lowland woodland plots. These plots were characterized by P. glandulosa, C. hookeri and Z. fagara and were concentrated on the higher pH, higher sand end of the axis. Plots from playa locations were quite distinct from them. Characterized by A. farnesiana. Playa plots occurred on the lower pH, higher clay, and higher SOC content end of Axis 1. The second CCA axis was primarily related to variation in soil BD (Figure 3A). Plots from grove and woodland communities exhibited the greatest range of variation along this axis (Figure 3B). C. pallida plants associated with soils having a relatively low BD and A. rigidula plants associated with soils with high BD defined the extremes of Axis II (Figure 3).

Figure 3

CCA Ordination (canonical correspondence analysis) results ordered in multivariate space along the first two canonical axes, separately depicting relationships among (A) tree-stature species (Table 1) and soil variables; and (B) communities in each 6 m × 6 m sample plot (symbols indicate community type) along a hill-slope transect.

Species centroids (+) indicate center of distribution among sample plots for each tree species coded as follows: acf = Acacia farnesiana; acr = Acacia rigidula; cel = Celtis pallida; con = Condalia hookeri; dio = Diospyros texana; kar = Karwinskia humboltiana; pro = Prosopis glandulosa; zan = Zanthoxylum fagara. Lines are vectors indicating direction of increasing value (from center outward) for soil variables as follows: BD = bulk density; Sand and clay = soil particle percentages; SOC = soil organic carbon; TN = total nitrogen.

CCA ordination of shrub-statured species exhibited greater overall variation (among more species) than tree species ordination. Total inertia was 7.85, of which 9.7% was explained by the first two axes (Axis I = 6.0%, Axis II = 3.7%). As with trees, Axis I was a synthetic gradient strongly associated with soil texture and pH (Figure 4A). The overlap among shrub cluster, grove and woodland plots observed for trees was also observed for shrubs; and, as with trees, plots from these communities were quite distinct from those in playas (Figure 4B). The centroids of most shrub species were relatively close to the plot origin, indicating a weak correlation with explanatory variables. Exceptions included P. aculeata, A. farnesiana and Z. obtusifolia in playas and some woodland plots (Figure 4). These species were positively associated with soils having higher clay content, and negatively related with soils having higher sand and pH. Axis II of the shrub ordination was mainly correlated with SOC, TN and BD (Figure 4A). Plots from woodland locations exhibited the greatest range of variation along this axis (Figure 4B); while plots from cluster and grove were more associate with higher soil BD, lower SOC and lower TN. Species like C. pallida, C. hookeri, and K. humboltiana were most strongly associated with soils having higher SOC and TN, whereas other shrub species (e.g., C. texensis, E. texana, A. greggeii) were more strongly associated with soils having relatively high BD.

Figure 4

CCA Ordination (canonical correspondence analysis) results ordered in multivariate space along the first two canonical axes, separately depicting relationships among (A) shrub-stature species and soil variables (Table 1); and (B) 2 m × 2 m sample plots (symbols indicate community type) along a hill-slope transect.

Species centroids (+) indicate center of distribution among sample plots for each shrub species coded as follows: acf = Acacia farnesiana; acg = A. greggeii; acr = A. rigidula; alo = Aloysia gratissima; bem = Bernardia myricaefolia; ber = Bernardia myricaefolia; cel = Celtis pallida; col = Colubrina texensis; con = Condalia hookeri; dio = Diospyros texana; eph = Ephedra antisyphilitica; eys = Eysenhardtia texana; for = Forestiera angustifolia; gym = Gymnosperma spp.; kar = Karwinskia humboltiana; par = Parkinsonia aculeata; pro = Prosopis glandulosa; sch = Schaefferia cuneifolia; zan = Zanthoxylum fagara; ziz = Ziziphus obtusifolia. Lines are vectors indicating direction of increasing value (from center outward) for soil variables coded as follows: BD = bulk density; Sand and clay = soil particle percentages; SOC = soil organic carbon; TN = total nitrogen.

Variation partitioning

Soil texture (% sand, % clay) accounted for the largest percentage of the unique total variation explained (TVE) for both tree (27.4%) and shrub (22.8%) species (Figure 5). A smaller percentage of the unique TVE was explained by SOC and TN (14.7% for trees and 18.7% for shrubs; Figure 5). However, the largest portion of the explained variation in both tree (57.9%) and shrub (58.5%) species distribution was that shared between both groups of soil variables (i.e., each group of variables redundantly explained the residual percentage of the TVE after the variation unique to each group was accounted for).

Figure 5

Variation decomposition shown as a percentage of total variation explained (TVE) for separate analyses of tree (Figure 3) and shrub (Figure 4) plots, that were uniquely related to variable groupings that included (a) soil organic carbon and total nitrogen, (b) % sand and % clay; or (c) an equal sharing by a and b.

Discussion

Isotopic evidence indicates the discrete cluster, grove and woodland communities at this site have all developed on what were grassland communities; and lines of evidence indicate these woodland communities are relatively recent (< 60-100 y of age)[1,46]. Thus, the discrete upland cluster, grove and drainage woodland communities represent examples of topoedaphic mediation of woody plant encroachment and community development. In this study we sought to determine if there were fundamental differences in the species and PFT composition of these communities that might shed light on how the shrub encroachment process plays out on upland and lowland portions of grassland landscapes. Our assessments of species and PFT composition and abundance were conducted across a local landscape-scale gradient. As such, they are not confounded by differences in climate, weather or land use history that can occur when comparing species abundance and community structure at multiple sites within a region or across bioclimatic zones. Here we only selected several categorical functional traits to define different PFTs, which we believe represent a wide range of functional traits that are important. However, including more detailed quantitative traits such as seed size or leaf traits may make our case stronger by providing other important aspects of the community assemblage.

Soil properties differed substantially between upland shrub cluster and grove communities [1] and along upland-to-lowland elevation gradients (Figure 1). Even so, upland cluster and grove communities were highly similar to each other and to lowland woodland communities (Figures 3a and 4a). Playa savanna communities, which are subject to periods of inundation that can cause woody plant mortality following heavy rains [34,35], were quite distinct from these, more so with the tree component than with the shrub component. However, with the exception of P. aculeata, the woody species and PFTs in playa communities were also present and typically common in cluster, grove, and woodland communities (Tables 2, 3, and 4). Furthermore, the extremes of the Axis 1 ordination gradient for both tree- and shrub-statured woody plants were characterized by similar PFTs (deciduous, N2-fixing, intermediate-to-deep rooting species; Figures 3, 4). Thus, we found little evidence to suggest that the major encroaching tree/shrub species or the PFTs they represent had strong preferences for or were uniquely confined to a given landscape location. Furthermore, it appears that patterns of both species presence/absence and their abundance along the hill-slope gradient were primarily a function of differences in intrinsic soil physical properties (e.g., texture and bulk density) and secondarily related to SOC and TN (Figure 5).

The woody communities on this former grassland site are relatively young (< 60-100 y of age) and may still be developing. Their species and PFT composition in the future may, to a large extent, depend upon the extent to which the present-day dominants are represented in the understory shrub layer. Overstory Prosopis plants (arborescent, deciduous, and N2-fixing) have been dying in cluster communities [1] as was evidenced by its absence in our sampling (Table 2). A. farnesiana, representing the same PFT as P. glandulosa, also occurred in clusters, but was not represented in the shrub layer (Table 3) and hence may not persist in the community. However other N2-fixing woody species occur in the shrub layer, so this PFT could still be represented in future vegetation states even with the demise of the current N-fixing PFTs in the tree layer. Prosopis overstory plants were doing well in the other communities (Table 2); and, with the exception of playa communities, was well-represented in the shrub class (Table 3). This suggests that overstory trees of this PFT will be replaced when they die. These observations centered around N2-fixing, deciduous PFTs are consistent with ordination results (Figures 3 and 4) and suggest two things: (i) significant changes in the species composition of each of these woody plant communities may be forthcoming; and (ii) species-specific rather than PFT-specific differences are driving the compositional dynamics of the woody communities developing on grassland. In other words, while some species may eventually be lost, the PFTs they represent will not be.

The N2-fixing PFT was a small component of the pool of encroaching woody plant species at our subtropical site from both a species richness and dominance (relative basal area) perspective (Table 4A and B). This is consistent with observations in Africa [47], despite the fact that leguminous trees are often regarded as a key component of African savannas [48]. However, a small contribution of N2-fixing shrubs to the species pool or total basal area, does not necessarily reflect their ecological importance. Indeed, early investigations at our study site suggested that non-N-fixing shrubs colonize soils that have been enriched in organic matter and nitrogen by N-fixing PFTs [49].

Rooting depth is an important functional trait in savanna ecosystems because niche separation of soil water use at different depth was believed to be an importance mechanism to tree-grass coexistence [49,50]. In our study site, as woody invasion is initiated by the establishment of P. glandulosa (deep rooted) on this former grassland, P. glandulosa can quickly grow roots deeper than the zone utilized by grasses, which enhanced their survival during the earlier stage of their life history [51]. Later in developed woody communities, hydraulic lift by P. glandulosa is an important mechanism of vertical water redistribution. However, the hydraulic lift is temporal dynamic and its effect on understory shrubs varies [31]. Overall, species interactions between overstory trees and understory shrubs are competitive and may cause the demise of many P. glandulosa trees [30,52]. Our observations of rooting depth PFTs in woody communities along the topoedaphic gradient (Table 4) suggest that the PFTs will very likely persist even though the dominance of individual species may vary in the future.

Relative basal area of evergreen and deciduous woody species are relatively stable along the topoedphic gradient (Table 4). Facultative evergreen functional type is more abundant in woodland and playa than the evergreen PFT, which is probably due to the associated soil moisture conditions along the edaphic gradient. It is very likely that deciduous and evergreen (including facultative) PFTs will persist in the communities with relatively stable percentages. However, more quantitative ecophysiological leaf traits (such as photosynthesis rate and specific leaf area) may be more indicative on how the functional traits may change over time, due to the differences between native and invasive species [21,53].

The number of woody species in shrub cluster communities and their patterns of occurrence at this site are strongly related to the size/age of the overstory trees [54]. Soils and microclimate change as the woody communities develop [1,55], and shrub species within these communities differ with respect to leaf longevity, specific leaf area, texture, and N content [21], daily and seasonal patterns of photosynthesis and water relations [30,32], functional rooting depths [31,52], and nitrogen responses [22,37]. These observations suggest a potential basis for differentiation among PFTs as woody plant communities develop on grassland. However, when confounding differences in abundance among shrub species are accounted for and compared against neutral predictions, random processes account for nearly all patterns of species and PFT occurrence [56]. This is consistent with observations that the classification of plants into functional groups based on leaf traits and N2-fixation potential are not necessarily reliable indicators of ecosystem processes related to resource use [21] or decomposition [57].

While there were numerous N2-fixing PFTs at our site, P. glandulosa was clearly the most aggressive initial invader of the historical grasslands [54]. This may reflect the fact that this species (i) is widely and effectively dispersed by grazing livestock [58]; and (ii) develops a root system that accesses water beyond the rooting zone of grasses very early in its life cycle [59]. Thus, traits related to seed dispersal and early taproot elongation may be more important in explaining the success of this species than N2-fixation. In addition, most of the other non N2-fixing woody PFTs at this site are bird-dispersed; and the establishment of livestock-dispersed P. glandulosa provides a perching structure attracting birds disseminating seeds of other woody PFTs [1]. This suggests that modification of soils by P. glandulosa plants colonizing grasslands may therefore be a secondary or coincidental factor influencing the encroachment of other woody PFTs – a conjecture consistent with our finding that patterns of species and PFT presence/absence and abundance were primarily a function of differences in intrinsic soil physical properties and secondarily related to modifications of SOC and TN by early-establishing shrubs (Figure 5).

Modification of soil properties beneath woody plant canopies is well-documented for this site [9,23,60] and for many other arid and semi-arid ecosystems undergoing increases in woody plant abundance (e.g., 61-63). Although several studies have implicated changes in soil C and N as a driver of plant distribution and species turnover [64-67], our results suggest that in the context of woody plant encroachment into grasslands, these effects are secondary to the influences of landform or geomorphic variables. However, in drier sites, shrub-soil feedbacks could play a more important role in driving woody encroachment. Variance partitioning results indicate that changes in SOC and soil N are likely to influence productivity and nutrient availability and may, over longer time-frames, feed back to influence PFT distributions. The two groups of variables also shared ~ 58% of TVE. This shared variation may reflect important interrelationships between relatively slow- and fast-changing edaphic properties or that both are related to another unmeasured variable. Both groups of variables exhibited similar spatial trends along the hill-slope, with higher values of SOC, TN, and clay in drainage woodlands than that in uplands (Figure 2). This suggests that correlations between these two groups probably contributed a significant portion of the shared variation. Dynamic simulations indicate that shrub-induced changes in SOC and TN will continue to increase for decades to come on this site [68]. However, while woody species composition may change, our data provide no basis for expecting that some woody PFTs will fare better than others.

The results from our study are consistent with assertions of random community assembly [56] in that woody species on this site were, for the most part, ubiquitously distributed along the topoedaphic gradient with its wide-ranging soil properties. Shrub species and the PFTs they represent in the Argentinian Caldenal, which is physiognomically similar to the Tamaulipan thornscrub of our Southern Great Plains study site, are also widely distributed across broad environmental gradients [69]. Thus, although the early establishment of Prosopis (a deep-rooted, N2-fixing, deciduous arborescent) appears to be important in initiating the transition from grassland to shrubland across the catena gradient at this site, our survey gives little reason to expect that the shrubland communities that subsequently develop are following any particular or predictable assembly rules that might consistently reflect PFT responses to changes in microclimate and soils occurring in the transition from grassland to shrubland or woodland. The lack of strong patterns in woody species or PFT occurrence within the four woody plant communities developing on contrasting topoedaphic settings may reflect the fact that PFTs can exhibit considerable convergence in resource use [70] and that there are many viable strategies for coping with given sets of environmental conditions. In the context of the physiognomic conversion of grasslands to woodlands, the ubiquitous distribution of shrub and arborescent species and PFTs on this site suggests that woody plants, and the diverse growth forms they include, are well-suited to a broad range of grassland topoedaphic settings.