Reconstruction of Piñon–Juniper Woodlands in the Sky Islands of the Davis Mountains, Texas, USA

: Piñon ( Pinus spp. L.)–juniper ( Juniperus spp. L.) woodlands’ historical stand structures were recreated to provide reference conditions and document long ‐ term changes in the Sky Islands of the Davis Mountains, Texas. Restoration of these isolated woodlands requires insights into the range of variability in current and historical stand structures, as well as an understanding of the spatiotemporal establishment and recruitment patterns of tree species. With drastic changes in forests and woodlands of the Southwestern United States widely reported, the main objective of this study was to reconstruct woodland tree temporal and spatial establishment patterns. A stratified random sampling approach was used to select two study sites each of 3600 m 2 in area. Within each site, all individual trees were mapped, measured, and cored for age determination. Age and tree location data were used to recreate the spatiotemporal patterns of tree species establishment and recruitment. Increments in density of both Mexican piñon ( Pinus cembroides var. cembroides Zucc.) and alligator juniper ( Juniperus deppeana var. deppeana Steud.) reached 422 trees ha − 1 in the 115 ‐ year period between 1890 and 2005; a yearly increment of 4 trees ha − 1 that did not reflect a rapid rate of change in these piñon–juniper woodlands. Age distributions reflected the multi ‐ cohort nature of these woodlands, and spatial autocorrelation measures were useful in the objective delineation of these cohorts. Temporal and functional niche differentiation of juniper was reflected in the development pattern where alligator juniper served as a pioneer species, exhibited a longer period of substantial recruitment, and had greater recruitment rates than that of Mexican piñon. Recruitment of Mexican piñon and alligator juniper occurred in an episodic fashion, with the majority of recruits being acquired between 1890 and 1949.


Introduction
Reconstructive studies enable the assessment of long-term changes and ecosystem dynamics, aid in the interpretation of ecosystem disturbance history, and are essential in guiding conservation and management goals [1,2]. Examination of changes in forest and woodland stand structure over time can be accomplished through the use of permanent plots that are monitored for tree recruitment, growth, and mortality at regular intervals. However, this approach may be limited in application due to inherent costs, as well as changes in land ownership. An alternative approach is to determine age structure retrospectively using dead and live trees [3]. The use of static age structures (at one point in time) in the interpretation of stand dynamics has been scrutinized, however. The absence of trees in past age classes has been posited to be a direct result of high rates of mortality and decomposition in lodgepole pine (Pinus contorta Dougl. ex Loud.)-Engelmann spruce (Picea engelmannii Parry ex Engelm.) stands, rather than lack of recruitment, as is commonly interpreted using static age distributions [4]. To circumvent the effect of a presumably high tree decomposition at early stages of stand development [5], it was recommended to construct age distributions at multiple points in time using fixed intervals [4]. This notion has been challenged, however, and static age distributions of ponderosa pine (Pinus ponderosa Lawson & C. Lawson) were deemed a faithful reflection of actual long-term establishment patterns rather than missing trees that died and decomposed with no trace, especially when standing dead and downed trees were included in these studies [6]. Evaluation of stand dynamics of a Mexican piñon pine (Pinus cembroides var. cembroides Zucc.) forest at the southern edge of its distributional range [7] found a preponderance of small trees (diameter <20 cm and height <6 m). The presence of two major cohorts and the lack of trees in past age classes suggested that recruitment and regeneration patterns were in discrete intervals rather than a continuous manner.
Sky islands are defined as an inland terrain that is made up of a widely spaced sequence of valleys and mountains [8], supporting a diversity of ecosystems and species within a few kilometers and relatively small changes in elevation. These isolated ecosystems, including the piñon-juniper (Pinus spp. L. Juniperus spp. L.) woodlands of the Madrean Archipelago, may be the least studied of all ecosystems in North America [9]. The Madrean Archipelago is characterized by high biological diversity due to the presence of two major floristic realms, the Rocky Mountains and the Sierra Madre Occidental, in addition to the convergence of three major climatic zones (temperate, subtropical, and tropical) [8]. The forest and woodland vegetation of the Madrean Archipelago has been classified [10] into three types: (1) encinal woodlands, (2) piñon-juniper woodlands and (3) high elevation coniferous forests.
There is a growing interest in restoring the forests and woodlands of the Southwestern United States [6,[11][12][13][14][15][16]. Restoration efforts will be an increasing challenge with climate change [17][18][19][20]. While Tausch et al. [20] and Miller et al. [21] reported an increase in tree density, accompanied by a decline in herbaceous production, as well as vertical expansion of piñon-juniper communities within the Great Basin, Blackburn and Tueller [22] noted the expansion of piñon and juniper into black sagebrush (Artemisia nova) communities in East-Central Nevada, and Landis and Bailey [15] reported substantial increments in density, basal area, and canopy cover of piñon-juniper communities since 1860. These changes in piñon-juniper communities were mainly attributed to overgrazing, fire exclusion, changing climate, and logging [21][22][23][24]. Restoration of these ecosystems requires the determination of reference conditions that can be used as approximations in guiding conservation plans and restoration treatments and take into account keystone variables capable of reflecting the evolutionary environment of a species or group of species [23,25]. These reference conditions can be defined using the range of historic variability in ecological structures and processes, which reflects the recent evolutionary history of these ecosystems [26]. Unfortunately, some restoration efforts were the product of a misconception in which the entire piñon-juniper type was considered "degraded" [27]. Although changes in structure and distribution have occurred throughout the piñon-juniper's extensive distribution range [15,20,22], some areas did not witness substantial increments in tree density or spatial extent [27]. For example, Ffolliott and Gottfried [28] reported an increase of approximately 148 trees ha −1 in 52 years (between 1938 and 1991) near Flagstaff, Arizona.
This study was conducted to provide reference conditions for restoration treatments and document long-term changes in stand structures of piñon-juniper woodlands of the Davis Mountains Sky Islands. The main objective of this study was to recreate historical stand structures by reconstructing tree temporal and spatial establishment patterns at five-points in time starting from the year 2005 to the earliest establishment date (i.e., 1985, 1965, 1945, 1925, and 1905).

Materials and Methods
The study was located within The Texas Nature Conservancy (TNC) Davis Mountains Preserve in Jeff Davis County, Texas. Climate is cool temperate with warm summers. Average annual rainfall is approximately 406 mm, with summer monsoons comprising two-thirds of the annual rainfall.
Mean summer temperature is 27 °C, whereas winter average temperatures range between −1 and 10 °C. Soils of the study area belong mainly to the Puerta-Madrone association (Alfic Lithic Argiustolls), with Loghouse association soils (Udic Haplustalfs), along stream channels of the study area [29].
The TNC Davis Mountains Preserve was stratified subjectively but without preconceived bias [30] into homogeneous areas based on vegetation, slope, elevation, and aspect. Stratification was achieved by merging 1 m Digital Orthophoto Quadrangles and 30 m Digital Elevation Models using a suite of tools in the ArcToolbox of ArcGIS 9.0 (ESRI, Redlands, CA, USA). To ensure that representative piñon-juniper communities were selected, areas less than 50 ha in size were excluded. From this sampling frame, two study sites (60 × 60 m) of similar vegetation and topography were randomly selected.
For each site, all overstory (height ≥1.4 m) vegetation was identified to the species level, and diameter at breast height (1.37 m; dbh) measured to the nearest 0.1 cm. All trees greater than 15.0 cm in dbh were cored at the base using an increment borer, and root collar diameter (RCD) and total height were measured to the nearest 0.1 cm and 0.1 m, respectively. Due to the high density of smalldiameter (dbh ≤15.0 cm) trees, a sub-sample consisting of 15 trees (1 tree per diameter class) per species was randomly selected for age determination. Within this sub-sample, trees >5.0 cm dbh were cored at the base, whereas trees ≤5.0 cm in dbh were cut and cross-sections collected. Dead, standing and downed, trees as well as stumps were included, regardless of decay stage, in the sample to extend age structure and minimize bias toward younger trees [4,6]. Dead wood accounted for 2% (8 trees ha −1 ) and 7% (67 trees ha −1 ) for site 1 and site 2, respectively. All trees within a site were positioned using a global positioning system receiver (Trimble GeoXT ® , Sunnyvale, CA, USA). For trees along the perimeter of the site (perimeter trees), distance to the nearest neighbor was measured to the nearest 0.1 m when the nearest neighbor was outside the perimeter of the plot [31]. For multistemmed trees, only the largest diameter (dbh) stem was considered for sampling. All seedlings (0.1 m ≤ height < 1.4 m) within each site were tallied and identified to the species level.
Tree cores were air dried, mounted, and sanded using progressively finer sand paper (50, 240, and 320 grit) following standard methodology [32]; cross-sections were processed in the same manner. Cores from the same species and same site were cross-dated (matching ring pattern) to ensure identification of missing and false rings by graphing ring-width measurements against the number of years [32]. For dead trees and stumps, death or cut dates were determined as the date of the outermost ring [33,34]. Decay was not sufficiently severe to prevent direct measurement and estimation of most death dates. Annual growth increments (ring-width) were measured to the nearest 0.001 mm using a Velmex TA system (Velmex Inc., Bloomfield, NY, USA). The quality of cross-dating was checked using COFECHA software [35,36].
Dbh and age frequency distributions were constructed and tree density (trees ha −1 ) and basal area (m 2 ha −1 ) calculated. Age structure (frequency distribution of tree age data) was constructed at multiple points in time using 20-year intervals. Estimates of tree density (trees ha −1 ) at multiple points in time were derived as the number of established trees that were still alive at that point in time.
To examine contemporary and historical tree spatial pattern, distance to the nearest neighbor for each tree within each site was calculated using an Arc Macro Language (AML) program in ArcInfo Workstation (ESRI, Redlands, CA, USA). Nearest neighbor distances were replaced with field measurements for perimeter trees, whose nearest neighbor lay outside the site. Nearest neighbor distances were pooled to calculate the observed mean nearest neighbor distance for each site [31]. The observed and expected mean nearest neighbor distance as well as the R statistic (the ratio of the observed mean nearest neighbor distance to the expected mean nearest neighbor distance) and the standard normal deviate (Z score) were calculated using Average Nearest Neighbor in ArcToolbox of ArcGIS 9.0 (ESRI, Redlands, CA, USA). The null hypothesis "the observed tree spatial pattern is random across the study site" was tested using an alpha level of 0.05 (−1.96 ≤ ZCritical ≤ 1.96).
For each site, a total of 20 adjacency matrices (adjacent tree pairs represented by 1 and nonadjacent tree pairs represented by 0) were constructed. Adjacency rule differed for each adjacency matrix and was defined in 3 m distance classes [37]. Adjacent tree pairs were defined as those pairs that are 0.1-3.0, 3.1-6.0, 6.1-9.0, and 57.1-60.0 m apart (one adjacency rule for each adjacency matrix). Adjacency matrices were used as the weighting matrices in calculating Moran's I spatial autocorrelation statistic for each distance class (d) for the variable tree age. Moran's I values were calculated using Spatial Autocorrelation in ArcToolbox of ArcGIS 9.0 (ESRI, Redlands, CA, USA). To test the null hypothesis "the observed spatial distribution of tree ages for each d is random", the standard normal deviate of that distance class (Z (d)) was compared to the critical standard normal deviate (−1.96 ≤ ZCritical ≤ 1.96) using an alpha level of 0.05. Correlograms were constructed by plotting the values of Z (d) on the vertical axis Y-axis) against d on the horizontal axis (X-axis) to explore the presence of even and uneven aged tree patches [37]. To identify even and uneven aged tree patches, Anselin Local Moran's I was calculated using Cluster and Outlier Analysis in ArcToolbox of ArcGIS 9.0 (ESRI, Redlands, CA, USA).

Results
Diameter distributions for all trees within each site were positively skewed, indicating an abundance (92%) of small-diameter (dbh < 20.0 cm) trees. Tree density and basal area were 711 trees/ha and 12 m 2 ha −1 , respectively. Alligator juniper was the most dominant species with mean basal area of 6 m 2 ha −1 and tree density of 283 trees/ha. Gray oak and Mexican piñon were ranked second in abundance with mean basal areas of 3 m 2 ha −1 and tree densities of about 200 trees/ha. Emory oak was a minor component but dominated the seedling strata.
Combined age distribution for piñon and juniper was unimodal, with a substantial establishment occurring from 1925 to 1945 ( Figure 1). Except for the period 1900-1929, recruitment rates of alligator juniper were always greater than piñon, with both species showing substantial increases in tree densities over time. Alligator juniper density increased from 33 trees/ha in the 1890s to 274 trees/ha in 2005, whereas Mexican piñon density increased from 7 trees/ha in the 1890s to 188 trees/ha in 2005 (Figure 1). Contemporary spatial patterns of Mexican piñon, alligator juniper, oak, and all trees combined were clustered and significantly deviated from a random distribution. Historical and contemporary piñon-juniper joint patterns were also clustered ( Table 1; Figures 2 and 3). Mexican piñon spatial patterns were significantly clustered, except in 1905, where the pattern was random (Table 1); alligator juniper spatial patterns were always significantly clustered.   Positive spatial autocorrelation for Mexican piñon tree age was shown for 2005, 1985, and 1965 for both sites ( Table 2). Positive autocorrelation indicated that tree clusters identified using nearest neighbor distances were of relatively similar age. Although not statistically significant, age spatial pattern had a negative autocorrelation for the periods 1925 for site 1 and 1905 for site 2, indicating that early establishment had trees with relatively contrasting ages during these times. Age data did not extend beyond 1905 for either site or species. Alligator juniper age spatial autocorrelation had a similar pattern to that of Mexican piñon with significant positive autocorrelation for 2005, 1985, and 1965 (Table 2). No significant spatial autocorrelation was revealed for alligator juniper age patterns at any time period for site 2. Correlograms, using adjacency matrices to define tree neighborhood, showed significant positive spatial autocorrelation at distances of 0.1-3.0, 6.

Discussion
Using diameter as a surrogate for age, reverse "J" shape diameter distributions are considered an indication of uneven-age, all-age, or multi-cohort stands [5]. Hump-like diameter distribution patterns may indicate the presence of at least two cohorts: an old cohort represented by large diameter trees and a young cohort represented by small-diameter trees.
Contemporary age structure indicated the presence of at least two cohorts: a remnant or initial alligator juniper cohort and a mixed piñon-juniper secondary cohort. The descriptors remnant and initial refer to two possible scenarios that may describe the development of this cohort. The initial scenario suggests the encroachment of alligator juniper into an open area, which may be attributed to comparatively greater seed longevity [38], greater tolerance of photosynthetic apparatus to water stress [39], greater drought tolerance, and better use of shallow soil moisture than piñon pine [40]. These pioneer members may have served as nurse trees for subsequent recruitment of piñon and juniper and may have also served as seed sources [38]. The remnant scenario, on the other hand, suggests opening of the stand by means of disturbance (e.g., firewood harvesting and/or fire) prior to 1890 that resulted in survival of this cohort. This scenario is consistent with the establishment of nearby Fort Davis in 1854. The secondary cohort corresponded to substantial increments in piñonjuniper recruitment over several decades . Substantial recruitment rates of Mexican piñon coincided with that of alligator juniper, although alligator juniper recruitment rates extended a little longer . The idea that junipers persist in the seed banks and contribute a continuous input into the seed bank supports this notion [38]. The period of substantial recruitment rates could be a result of favorable establishment and survival conditions and/or a result of response to disturbance in the form of fire exclusion and/or overgrazing [20,22,24] [21,23,24]. A portion (1935)(1936)(1937)(1938)(1939)(1940)(1941)(1942)(1943)(1944)(1945)(1946)(1947)(1948)(1949) of the period of substantial recruitment rates for piñon-juniper exhibited an above average annual precipitation, and the period was preceded by the establishment of Fort Davis and the Euro-American settlement of the area. This establishment period coincides with establishment periods for piñonjuniper in Northern Arizona [15] and the Great Basin [20].
Piñon-juniper age distributions between 1905 and 1945 suggest a stand initiation stage which occurred over several decades. The distribution of Mexican piñon indicated less variation in age as compared to alligator juniper, and this may be attributed to shorter recruitment periods of Mexican piñon. In addition, age distributions suggested a longer lifespan for alligator juniper, which may also have served as a pioneer species in structural development.
Greater recruitment rates and density estimates of alligator juniper between 1890 and 2005 could be explained by the same factors potentially responsible for the establishment of junipers in canopy interspaces and open areas [38][39][40].
Significant spatial clustering may be explained by seed dispersal mechanisms in which piñons are mainly dispersed by birds (e.g., piñon jays (Gymnorhinus cyanocephalus)) and rodents (e.g., piñon mice (Peromyscus truei)),whereas junipers are mainly dispersed by birds (e.g., bluebirds (Sialia mexicana and S. currucoides)), rodents, as well as mammals (e.g., mule deer (Odocoileus hemionus) and coyote (Canis latrans)), after passing through the digestive system [38]. The spatial clustering may also have created habitats for this variety of animals, further enhancing seed dispersal by these species. Animal-dependent seed dispersal (scattering in multi-seed caches), seed dispersal by gravitational forces, and the spatial heterogeneity of the environment may have contributed to the observed spatial pattern.
The observed spatial pattern may be confounded by the inability to capture spatial arrangement of seedlings that might have occurred in great numbers and experienced high rates of mortality and decay [5]. Temporal spatial patterns indicated that Mexican piñon and alligator juniper were equally successful in inter-canopy and open area establishment as well as establishment within existent clusters, in contrast to Everett et al. [41] who found that piñon seedlings were more abundant under tree canopies than in inter-canopy spaces in Utah. In the Davis Mountains, however, summer monsoons may have lessened the need for nurse-tree micro-environments that are necessary for establishment and survival in the Great Basin [38]. In addition, differences among species (Mexican piñon and alligator juniper versus piñon pine and Utah juniper) in terms of light, moisture, competition tolerance, and nutrient levels needed for germination and seedling survival and development may have attributed to the observed differences. Climate change would certainly impact these development patterns, and under future climate projections, nurse-tree microenvironments may become more important for establishment and survival of piñon and juniper in the Davis Mountains [17][18][19].
Significant positive spatial autocorrelation of Mexican piñon in 2005, 1985, and 1965 indicated the presence of piñon clusters of relatively similar age; this similarity may be attributed to the lack of piñon density variation among these periods. In 1945In , 1925In , and 1905, the spatial distribution of ages was similar to that expected from an independent random process. The presence of clusters of relatively similar age piñons and clusters of relatively dissimilar age piñons reflected the multi-cohort nature of Mexican piñons.
Under the assumption that juniper regeneration occurred in a clustered pattern, differences in the spatial distribution of successful recruitment between 1925 and 1945 may be attributed to higher rates of seedling mortality (e.g., low intensity surface fire or a more pronounced within-cluster competition) prior to 1925. In general, the shift from negative to positive spatial autocorrelation emphasized the clustered recruitment patterns of piñon and juniper, and contemporary spatial distributions of piñon-juniper ages reflected the irregular (different size clusters) multi-cohort nature of these stands.

Conclusions
Stand structure reconstruction in this piñon-juniper woodland revealed substantial increments in density, similar to those reported in the Great Basin and Northern Arizona [15,20]. However, when viewed in the context of the elapsed time, annual increments were similar to those reported by Ffolliott and Gottfried [28] and did not reflect a rapid rate of change in these piñon-juniper communities. Age distributions reflected the multi-cohort nature of these communities. Temporal and functional niche differentiation of alligator juniper was reflected in the development pattern where alligator juniper served as a pioneer species, exhibited a longer period of substantial recruitment, and had greater recruitment rates than that of Mexican piñon, which may have been influenced by climate conditions that favored germination and establishment of woody species in this region. Recruitment of Mexican piñon and alligator juniper occurred in an episodic fashion with the majority of recruits being acquired between 1890 and 1949. Examination of spatial pattern and spatial autocorrelation measures suggested that recruitment occurred in different-size clusters within the canopies of existent clusters as well as among inter-canopy spaces.
Restoration efforts should consider other lines of evidence (e.g., historical photographs and early historical accounts), in conjunction with stand reconstruction data, when establishing reference or desired conditions for restoration treatments. Although stand structure reconstruction may have provided insight into the development pattern of these piñon-juniper woodlands, an understanding of the development pattern is far from being complete without consideration of its contextual past climate and disturbance history [42]. When stand structure reconstruction is coupled with the evaluation of disturbance history, especially the reconstruction of fire frequency, and tree species distribution climate change projections a better understanding of the development patterns and restoration of resilient woodlands would be attained.