Abstract

:
To evaluate the impacts of urbanization on soil and vegetation in protected forest areas, 12 forest sites in Southeastern Michigan USA were studied in an indirect urban–rural gradient. Field study plots were established in forest edge zones of each protected area. Significant findings were that in these edge zones of protected areas: (a) soil nitrogen tended to be greater where surrounding housing density was greater; (b) overstory woody biomass and basal area were greater where surrounding housing density was greater; and (c) the study region overall exhibited low soil carbon content (mean 2.71%) and relatively high soil nitrogen content (mean 0.20%), yielding a surprisingly low surface soil C/N ratio (mean 13.4). Overall, 24 woody plant genera were encountered, with the three genera Acer, Carya and Quercus accounting for 83.7% of total biomass and 74.1% of total basal area. No significant relationships were observed between housing density and soil C/N ratio or between housing density and foliar N. Results indicate that a halo of urban-ecological impacts exists in the landscape of Southeastern Michigan, similar to previously studied linear urban–rural gradients in other regions.

1. Introduction

High levels of human activity are commonly associated with ecosystem change [1,2,3,4], which makes the creation of nature preserves an important element of conservation planning within a human dominated landscape [5,6]. However, there is a growing appreciation that the ecosystem structure and function within a protected area can be affected by increasing urbanization influences in the surrounding landscape [7,8,9,10]. As a result of the shared hydrologic, biogeochemical, and atmospheric resources connecting protected areas and the surrounding urban landscape [11], traditional protection efforts may be insufficient to maintain ecosystem integrity within protected patches of land. Such efforts typically change direct human activity but do not directly address slow, long term changes to shared resources. As a consequence, even after setting aside land, urbanization in surrounding landscapes may interfere with conservation goals [12,13,14].

Conversion to residential land use, a particular alteration of the landscape, has been tied to specific changes in ecosystem processes [15,16]. At a minimum, increased management concurrent with urbanization can dictate vegetative biomass and N concentration across the landscape [17]. At higher housing densities, urban settings have also been linked to secondary impacts that include increased nighttime temperature, precipitation, and nutrient deposition [18,19,20]. Such urban-derived effects extend into nearby protected areas and alter ecosystem function [2,21,22]. As a result, urban presence has been identified as an important context for N cycling with potentially great impacts on urban vegetative composition [23].

Ecosystem preservation techniques suitable for an urbanizing landscape are of increasing importance. In the United States, the population is progressively moving away from metropolitan centers and into low-density suburban and exurban developments [24]. In the year 1950, only 6% of area in the conterminous United States was developed at a density of forty acres (16.2 hectares) or less per home [25]. By the year 2000, census block groups with densities of less than forty acres per home constituted 27% of the conterminous USA [25]. This change was primarily driven by expanding low-density development (1–40 acres/house) rather than intensifying urban areas (less than 1 acre/house (0.4 hectares/house)). Two consequences of the growth of low-density development since 1950 have been an increase in the number of natural areas with proximal human impact [26,27,28], and a related increase in the importance of understanding this impact. Despite documented relationships between changes in surrounding land-use and altered ecosystem function, many questions remain. Of particular relevance is the potential for urbanizing landscapes to affect the biogeochemistry and ecosystem function in remnants of wildland soil and vegetation [7,8].

Given the wide range of potential human impacts on protected ecosystems, Russell [28,29] divided the potential influences of human activity upon protected areas into the categories of “obvious” and “subtle” effects. “Obvious” effects constitute those human activities that have clear before and after components. By contrast, “subtle” effects normally operate at longer time scales relative to typical observation periods. Such effects can be more difficult to distinguish than the typically studied obvious effects. Consequently, subtle effects are commonly accepted as status quo rather than as anthropogenically driven changes.

Many impacts of urbanization on biogeochemical cycles fall into the category of “subtle” effects. Residential land-use may produce subtle effects such as long-term alteration of ecosystem function and compounding time-lagged effects [30,31]. Further, resultant biogeochemical changes are difficult to directly mitigate because of their scale, pervasiveness, and lack of controls on dominant drivers [32]. As an example, the construction of a local road network leads to generally permanent increases in N deposition [33,34]. Introduction of nitrogen (N) from the surrounding landscape to protected areas is a serious concern [35,36], as N is frequently the limiting element in terrestrial ecosystems [31], and N availability is often a determinant of species richness, diversity, and productivity [37,38]. The potentially significant effects of nitrogen addition on unmanaged ecosystems, along with increased human presence near protected areas, have fueled a number of long-term ecological studies on the effects of urbanization on biogeochemical cycles [24,37,39,40]. A recent biogeochemical field study in residential yards found that homeowner practices resulted in large fluxes of N to soils in residential neighborhoods [41], however, a recent chronosequence study of residential yards following conversion from other land uses concluded that residential soils can retain large quantities of C and N over time [42], raising the question of whether N additions in residential land would affect N cycling, and resulting changes in C storage, in protected areas nearby.

There is a longstanding interest in understanding relationships between land use/land cover and the storage of carbon (C) in soil and vegetation. Protected areas, by allowing large trees to grow and by minimizing soil disturbance, may be important locations for C storage in the landscape. In the urbanizing Front Range of Colorado, established urban green spaces were found to harbor larger C pools than native grasslands or agricultural fields [18]. New York City alone is estimated to store 1.2 million metric tons of carbon in soil and vegetation [43]. Similarly, In the Phoenix metropolitan region, human activities may have increased soil organic matter by more than 40% [20]. The few studies conducted to quantify C pools in urbanizing landscapes are unusual because ecologists and biogeochemists have historically focused more attention on wild land ecosystems. Little is known about the general patterns of carbon storage in urbanizing landscapes.

Here we studied the relationships of the response variables of soil and vegetation nitrogen pools, carbon pools, and vegetation biodiversity in protected areas to the potential driving variable of housing density in the surrounding landscape. We focused on protected forest areas and within those areas the edge zones adjacent to non-forest land cover. We hypothesized that soil and vegetation would be more N-rich in edges of protected areas surrounded by greater housing density, potentially as a result of urban pollution and runoff [11]. We also hypothesized that differences in N-richness would correlate with differences in vegetation biodiversity and pool sizes of C in biomass. In addition, by taking place in Southeastern Michigan, our study extends the questions and approaches of urban-rural gradient studies to a new region, the Upper Midwest United States.

Southeastern Michigan is an area well suited to address these questions because it has been undergoing development driven changes in land use/land cover for several decades [26] and it is subject to regionally elevated nitrogen deposition [44]. Existing long-term urban–rural research sites (Maryland, New York/Connecticut, and Arizona) are 1000–3000 km from the present region and have differing climates and soils. Southeastern Michigan’s spatial pattern of human housing density in the landscape also differentiates this study from previous studies. The study area (Figure 1) features three distinct population centers (Detroit, Jackson, and Ann Arbor) with spatially complex population gradients between urban and rural areas [45,46]. In typical studies along urban–rural gradients, population density follows a monotonic decline with distance from an urban center [47]. Here, with three spatially distributed urban centers and the lack of such a monotonic pattern, we generalize the concept of an urban-rural gradient to an indirect gradient. We use housing density as the independent variable, making our approach similar to previous studies along urban-rural gradients. This approach allows us to explore the application of the gradient concept in a landscape that contains a range of sizes of urban centers together with non-monotonic transitions in population densities between and around them. Southeastern Michigan has a range of housing densities spanning more than three orders of magnitude (Figure 1).

Figure 1.
Map of Housing Density, Major Cities, and Site Locations. Housing density in census block group resolution is shown, with the three major metropolitan areas and site locations used in the present study. Acres/Unit refers to average acres per US Census Bureau housing unit within each census block (1 acre/unit = 0.4 hectare/unit). Shapes of site symbols represent housing density classification used in the present study [48,49].

Figure 1.
Map of Housing Density, Major Cities, and Site Locations. Housing density in census block group resolution is shown, with the three major metropolitan areas and site locations used in the present study. Acres/Unit refers to average acres per US Census Bureau housing unit within each census block (1 acre/unit = 0.4 hectare/unit). Shapes of site symbols represent housing density classification used in the present study [48,49].

2. Experimental Section

Six contiguous counties totaling 10,522 km2 (Jackson, Livingston, Macomb, Oakland, Washtenaw, and Wayne) in Southeastern Michigan, USA composed the study area. This area contains the urban centers of Detroit (population 918,849), Ann Arbor (population 114,024) and Jackson (population 36,316), as well as large expanses of medium and low density development [48]. Since our focus was on forest patches, we identified protected areas containing forest area greater than one hectare [50]. Block groups from the 2000 Census [48] were used to stratify the study area to ensure that sites were selected across a range of housing densities (Figure 1). Stratification was accomplished through a proxy variable that represented land acres per housing unit (LA/HU, [51]), or the inverse of housing density. We contacted site managers about management practices and usage trends so as to avoid “obvious” anthropogenic impacts that might disqualify a site. Examples of disqualifying management activity included controlled burns and removal of invasive species in the area surrounding the plot. Site selection criteria included level of recreational usage, described management activities, and an informal site history. Most sites managers indicated that the forests were established on former agricultural fields and none contained old growth forest. Twelve sites, equally distributed between urban, suburban, exurban and rural categories of housing density (Table 1), were chosen so as to create an indirect urban–rural gradient throughout the six-county area [22].

All preserves met the minimum requirement of a closed canopy and touching crown branches of adjacent trees. The preserves were designated between 1903 and 2006, a difference of more than one hundred years. Within the 1000 m distance surrounding each preserve, the oldest housing development was completed between 1918 and 1970 (Table 1). Field work was conducted between June 1 and July 27, 2007.

2.1. Sampling Design

At each site, a single 20 m × 50 m (0.1 ha) research plot was established to quantify ground cover, identify woody plants, sample soils, and sample foliage. Plots were positioned within the forest patch with the outermost 50 m plot edge located along the boundary with non-forest land cover and the innermost 50 m plot edge located 20 m into the forest patch. We thus refer to these plot locations as forest edge zone. Edge zones were chosen for study because of their potential importance in retaining N in runoff from the surrounding landscape and because of their likelihood to reflect other impacts from the surrounding landscape such as invasive plant species.

Based on Bonham [52], ground-cover was defined as the percentage of ground surface covered by living plant material between the ground and a plane 1m in height. As ground-cover can be interpreted differently depending upon measurement technique [53], three different methods were used for recording ground-cover. Ten 5 m × 1 cm transects were set up along the center axis of the plot, ten 1 m2 (0.5 m × 2 m) subplots were set up along the plot borders, and two 10 m2 (5 m × 2 m) subplots were set in opposite plot corners. Using methods described below, presence and absence were noted for grass, sedge, woody plant, forbs, and obstructed cover [52]. Obstructed cover included stones or other non-living materials.

Ground cover in each of the 1m2 subplots was observed using a frame with border markings delineating 1% squares. A hypothetical 1 m column extending upwards from ground level in each 1% square was evaluated for presence or absence of each cover type. Each cover was separately quantified, allowing overlap, which resulted in an aggregate percent cover that frequently exceeded 100%. A similar method was used for transect measurements.

Once ground-cover measures were completed, all trees in the research plot with a diameter at breast height (DBH) greater than 1cm were measured for DBH and identified to genus. Using genus rather than species enabled rapid field assessment of woody plant diversity [54,55]. Five to ten sun-exposed leaves were then sampled from each of approximately twenty trees. Samples were clipped using pruning shears attached to a two meter pole. Efforts were made to ensure that leaf sample composition roughly reflected site composition, but genera such as Juglans and Quercus were occasionally underrepresented if there were limited leaves within two meters of the ground.

Four stratified randomly located, soil samples were taken at each site. Woody debris and leaf litter was cleared from the forest floor, and a 3 cm × 17 cm (diameter × depth) soil core sample was removed from the surface soil, containing the Oea horizon (soil layers with a high percentage of organic matter) where present and upper mineral soil. Each core was separately homogenized to mix horizons within the individual core. A similar method, including removal of loose leaf litter, was used in Michigan by Pregitzer et al. [56]. Soil samples were transported to the laboratory in airtight plastic containers.

2.2. Analysis of Field Samples

Immediately following collection, leaf samples were dried overnight in a forced air oven at 65 °C, while the homogenized soil was dried overnight in a gravity convection oven at 105 °C. Dried samples were stored in a climate controlled room for up to 3 months. Soil samples were then sieved (2 mm). Sieved soil samples were analyzed for texture using Bouyoucos hydrometer methods [57]. Soil and leaf samples were analyzed for carbon and nitrogen content by dry combustion in a Carlo Erba NC 2500 elemental analyzer attached to a Thermoquest-Finnigan mass spectrometer. Overall, 241 unique leaf samples and 48 soil samples were analyzed.

Stem densities of plants within the plot were calculated from field data. For woody plants, size, above-ground biomass (dry weight), foliar biomass (dry weight), and basal area were calculated from DBH using published allometric relationships by genus [58,59,60,61,62,63,64,65,66]. To evaluate compositional diversity, the Shannon Diversity Index was calculated from genus-level data at each site. While using genus-level data leads to lower measured diversity than species-level data, the genus is more easily assessed rapidly in the field and relative diversity estimations among sites are still useful [67].

Total foliar carbon and nitrogen at each site was calculated using foliar C and N concentrations and allometric equations. In order to use species-specific allometric equations with genus-level data, US Forest Service Forest Inventory and Analysis Program data on stem count were acquired for the study region [68]. The allometric equations did not uniformly allow separation between sun and shade leaves, so no such separation was made in scaling up foliar C and N. The relative proportions of the three most prevalent species of each genus in the six-county study area were used. The foliar biomass of each field-measured tree was calculated with each allometric equation. A composite foliar biomass was then created for each tree by proportionally weighting the species specific foliar weights by the species’ prevalence. Foliar mass was then summed by genus at each site. By site, foliar C and N concentrations were multiplied by the genus-level summed foliar mass to produce site totals of foliar C and N. Allometric whole tree biomass was calculated similarly, using whole tree biomass equations, and similarly summed to produce site-level tree biomass.

For each field site, an area-weighted average of surrounding housing density was calculated among the census blocks groups within a circle of a 1000 m radius from the exact center (as measured by GPS) of each field research plot. Using 30 m resolution digital elevation models (DEMs) site elevation, aspect and topographic index of wetness were calculated [69] using ArcGIS software (ESRI, Inc.).

2.3. Statistical Methods

The categorical classification scheme used for housing density (Figure 1) has categories with progressively wider ranges of density. With an equivalent number of research sites in each density category, the scheme produced a site distribution which was skewed towards the densest urban areas. To normalize the site distribution, a loge transformation was used on the housing density variable. This transformation improved the fit between the independent and dependent variables’ distributions and reduced rural sites’ leverage. Future references to “housing density” are equivalent to loge (LA/HU).

Moran’s I was used to test for any unaccounted factors causing spatial autocorrelation [70,71]. Tests of Moran’s I were run in STIS 1.62 [72] on mean %N in soil, C/N ratio in soil, %N in foliage, and C/N ratio in foliage. The oldest documented housing development age within a 1000 m distance from each field site was tested against housing density in a one-way analysis of variance. All remaining statistical analyses were run on SPSS 15.0 [73]. Summary linear regressions were used with measures taken once at each site. In order to incorporate the internal variability in measures that were taken repeatedly at a single site (such as soil N concentration), an additional multi-tiered analysis was pursued through a linear mixed model. The mixed model is useful for reducing type I errors in clustered data with potential causal heterogeneity [74,75]. The linear mixed model was run with heterogeneous compound symmetry covariance, which in this dataset implies that the relationship between sites varies, but within each site the samples share the same variance structure. Akaike’s Information Criterion (AIC), rather than R2, is typically used to compare relative performance of linear mixed models [76]. Significant differences in linear mixed model performance were measured through chi-squared goodness of fit tests.

3. Results and Discussion

3.1. Relationship between Soil N and Surrounding Housing Density

Soil N concentrations in forest edge zones across the 12 protected areas exhibited no spatial clustering or simple directional gradient across the study region (Moran’s I = −0.16, p = 0.16), but did demonstrate a weak and marginally significant relationship with housing density in the surrounding landscape (univariate linear regression, R2 = 0.27, p = 0.08; Figure 2). A linear mixed model was run to reduce the possibility of a type I error. Results from the linear mixed model indicate that incorporating housing density significantly improved the prediction of soil nitrogen concentration (χ2 goodness of fit test, p < 0.01) over the null model (Table 2 and Table 3). There was no significant relationship between the oldest housing development surrounding each site and housing density surrounding each site (F(3,8) = 0.79, p = 0.53) and thus housing age was not included in further analyses. In combination, the two statistical models provide evidence of a significant, albeit weak, correlation between soil N concentrations and surrounding housing density (Table 2 and Table 3).

Figure 2.
Soil % N. Mean soil N concentration in the 12 study sites versus mean loge land acres per housing unit (LA/HU) in the surrounding landscape. Note that higher values of housing density occur at the left of the diagram, at lower values of (LA/HU). 1 land acre per housing unit = 0.4 ha per housing unit. (a) Where surrounding housing densities were greater, mean %N in soils of protected areas tended to be greater (R2 = 0.27); (b) Site elevation was also a strong predictor of soil %N (R2 = 0.52, p < 0.01) across the twelve study sites.

Note: The relationship between housing density in the surrounding landscape and soil nitrogen concentration in a protected area is shown using a linear regression on the site mean values. Data were derived from 12 field study plots in the present study. Loge land acres/housing unit (LA/HU) is used as a proxy for housing density. Model coefficients, R2 and p values are listed for linear regressions.

The relationship between housing density in the surrounding landscape and soil nitrogen concentration in a protected area is shown using a linear mixed model on all samples. Data were derived from 12 field study plots in the present study. Loge land acres/housing unit (LA/HU) is used as a proxy for housing density. Model coefficients and p values are listed for these statistical tests. Akaike’s Information Criteria (AIC) allows relative comparison of linear mixed models. The χ2 goodness of fit test indicates whether a model is a significant improvement over another model. In this case, all models are compared against the null model.

3.3. Vegetative Composition of Sites

Of the 2027 trees that were identified across the study region (Figure 3), 24 total genera were encountered but the three genera Acer, Carya and Quercus accounted for 83.7% of total biomass and 74.1% of total basal area. Stem counts were dominated by Prunus, Quercus, Ostrya, and Acer, which together accounted for 67.3% of total woody stems. In the eleven genera with the least abundance, fewer than 25 stems (1% of total) were encountered. Shannon index values (mean of 1.52) among the sites were not statistically different, indicating relatively uniform biotic diversity in the woody plant communities across sites.

Figure 3.
Aggregate Overstory Composition Results: 2,027 woody plant stems (trees and shrubs) were encountered in the 12 study plots established in forest edge zones. One plot was established in each of the 12 protected areas throughout the study region. (A) Frequency distribution of measures of diameter at breast height (DBH); (B) Percentage of woody basal area by genus across all woody plants; (C) Percentage of woody biomass by genus across all woody plants; (D) Percentage of woody stem counts by genus across all woody plants.

Figure 3.
Aggregate Overstory Composition Results: 2,027 woody plant stems (trees and shrubs) were encountered in the 12 study plots established in forest edge zones. One plot was established in each of the 12 protected areas throughout the study region. (A) Frequency distribution of measures of diameter at breast height (DBH); (B) Percentage of woody basal area by genus across all woody plants; (C) Percentage of woody biomass by genus across all woody plants; (D) Percentage of woody stem counts by genus across all woody plants.

When the site with the fewest stems and lowest biomass (Black Pond Woods, BPW) is excluded, protected areas surrounded by lower density housing had significantly lower total overstory biomass (Figure 4) (R2 = 0.60, p < 0.01) and basal area (R2 = 0.44, p = 0.03). A hat matrix of this regression, to identify unduly influential or exceptional data points, indicates that the BPW site exhibited undue leverage on the regression results (Cook’s distance of 0.67, studentized residual of −4.40, [77]). Contrary to our hypothesis, the overstory biomass was not driven by soil N concentration (R2 = 0.02, p = 0.72), suggesting that the vegetative growth in the region may not be N-limited.

Relationships between overstory and understory vegetation have been seldom studied in human-dominated landscapes, where forest patches are fragmented and increased edge habitat could be expected to alter light availability and temperature, while providing an inroad for invasive plant species. We found that higher overstory tree biomass correlated with a higher proportion of bare ground cover (R2 = 0.67, p < 0.01) and a decrease in woody plant ground-cover (R2 = 0.37, p = 0.04), indicating that greater overstory biomass was associated with more open, less dense understories even in the edge areas that we studied (Figure 5).

3.5. Soil N Content and Housing Density in the Surrounding Landscape

Previous research in New York, North Carolina, and Arizona has suggested relationships between housing density and soil N, which are supported by the present findings of a weak relationship between housing density and soil N content in forest edge zones (Table 4). We found generally stronger relationships between soil N and physiographic variables, including elevation (Table 2 and Table 3). Over our range in elevation of about 350 m throughout the region, our finding of lower soil %N at higher elevations could be result of mobile N being carried by soil water toward streams and rivers. The low soil C/N ratio variance (= 1.5) at our sites suggests a high degree of uniformity across the region, and makes statistically significant relationships from anthropogenic influences difficult to detect in combination with other spatial variables. Even though statistically significant changes are not detected, biologically meaningful changes may be occurring.

This concern about marginally significant results is echoed in similar previous research, wherein results were discounted due to confounding factors. Examples abound: along a New York urban–rural gradient, a positive correlation between population density and N mineralization in forest stands was discounted due to earthworm activity [38,77,78]; a positive correlation between net N mineralization and urban land-use in North Carolina was criticized for insufficiently accounting for urban heat island effects [40,79]; and some research conducted in Arizona was unable to fully substantiate differences in N mineralization rates between urban and rural areas [80]. Taken together with the present study, these previous results suggest either there is no such effect, or that effects are subtle and difficult to disentangle amidst confounding variables.

Our results also showed that in Southeastern Michigan physiographic factors such as topography and soil texture are currently as important, or more important, in explaining landscape-scale patterns of soil N as housing densities. The present analysis is confounded by an apparent relationship between housing densities and topographic and soil factors, such that higher housing densities are more likely at lower elevations and with relatively lower sand content in soils of the region. These challenges are not unique to this study; complex entangled causes are endemic to much landscape level ecological research. However, a general preference for modeling biogeochemistry with soil and topographic inputs as opposed to causes related to human activity is consistent with other research in the field [78,40]. Continued research into the complex relationships between anthropogenic (housing density and the impacts that accompany human presence) and physiographic or topographic drivers is likely to be key in understanding the biogeochemistry of the human-dominated landscape. Another factor that could be important in affecting surface soil N concentrations is the distribution of vegetation species from one location to another. Different tree species could result from different site histories, ages, or planting, and could affect surface soil N concentrations through litter inputs.

Our method of mixing Oea and upper mineral horizons during sampling was designed to allow us to sample a large number of sites with limited intensity, as required in adapting methods from intensive studies to use in a landscape-scale study. This makes it difficult to directly compare our soil C/N ratios against those from intensive-study forest sites, as intensive studies generally collect separate samples of O horizons and mineral soil horizons. However, the unexpectedly and uniformly low site C/N ratios we found (mean = 13.4) make this a compelling post hoc question: are forest edge zone soils uniformly low in C/N ratios in our region? Additionally, open questions include whether C and N biogeochemistry in interiors of protected forest patches differ from that of forest edges, and under what conditions surrounding housing densities could show effects that penetrate to forest interiors. The present study focused on the forest edge in order to maximize the likelihood of detecting impact. Follow-up work could be conducted in interior areas of forest reserves, and to determine whether vegetation and soil C and N differ between the forest edges and interiors in these protected areas.

Biogeochemical and overstory summary data are presented for each of the 12 field study plots in the present study. Loge land acres/housing unit (LA/HU) is used as a proxy for housing density; lower values of (LA/HU) indicate greater housing densities (1 land acre per housing unit = 0.4 ha per housing unit). C and N respectively stand for carbon and nitrogen. Biomass is expressed as kg dry weight; biomass and basal area refer to all woody vegetation in one field study plot (0.1 ha) from each study site. Site abbreviations are based on park and town names.

The C/N ratio of soils in this study are comparable, but on the low end, of those measured in mineral soils of other forests in the N-rich northeastern United States [81]. A mean C/N ratio of 18.5 was found in the upper mineral soil in the White Mountain Region of New Hampshire [82], and Northern Connecticut forests possessed mean C/N ratios in upper mineral soil horizons from 15.5 to 18.9 [83]. McFarland [84] took soil samples of 20 cm depth and found C/N ratios ranging from 13 to 21 in several northern temperate forest ecosystems. In soil cores from four Northwestern Michigan Forests starting at 10 cm in depth and continuing to 30 cm in depth, Pregitzer [56] found C/N ratios ranging from 13.8 to 16.2. These soil samples are from a nearby region, exclude the Oea horizon, and still contain higher mean C/N ratios than found in the present work.

The types of sites investigated in these earlier studies, which were generally conducted in wildland forested landscapes, stand in contrast to our study of forest edge zones in a human-dominated landscape. Even so, these studies demonstrate the large-scale regional consistency of relatively low C/N ratios in upper mineral soils of the Northeastern US. Unfortunately, our measures of soil C/N are not directly comparable to these intensive forest studies because of differences in sampling methods, but indicate that further study is warranted to determine how N levels or C/N ratios in the soils in our region compare with N levels or C/N ratios in other areas of the Northeastern US. Several factors may drive low C/N ratios such as elevated soil N content, low soil C content, or the known presence of mineral soils in the samples collected for this study [83].

There is evidence that the foliage in overstory woody plants covered by this study showed moderately high levels of foliar N in comparison to other studies. Our mean N concentration across all overstory foliage = 2.14% (Table 4). For comparison, we analyzed side by side data for foliar N and net N mineralization (a key measure of soil N cycling rate) across 30 forest plots widely distributed across New Hampshire ([85], their Table 1). The 15 sites in the New Hampshire study with lowest rates of net N mineralization averaged 53.7 kg N ha−1 y−1 and had an average foliar N of 1.53%. The 15 sites in the New Hampshire study with the highest rates of N mineralization averaged ca. 3 times higher N mineralization, 118 kg N ha−1 y−1, and had an average foliar N that was also higher, at 2.09%. The median value of foliar %N across the New Hampshire study was 1.87% [85]. A study in the nitrogen-rich Adirondack Mountains (New York, USA) found mean foliar N concentrations that ranged from 2.50% to 2.96% [86]. The current study found foliar N concentrations (Table 5) above the median reported by Ollinger et al. [85] in New Hampshire but below those reported in the Adirondacks. Moderately high foliar N across the wide range of species we observed in Southeastern Michigan may be attributable to a number of factors. These include high rates of N deposition in the region [87] and the agricultural history of the region, among others.

Foliar sample data are presented with number of samples, mean % nitrogen by mass (ash-included basis) and mean C/N ratio by mass. Data were derived from 12 field study plots in the present study.

3.6. The Effects of Time and Agricultural Abandonment

Present day forests in our region are frequently located on former agricultural fields, which may be a contributor to the low C/N ratios [88]. At time of abandonment, an agricultural field typically has reduced (≈30%) carbon content as compared to uncultivated areas [89,90]. With carbon content a major factor contributing to nitrogen retention in soil [91], date of agricultural abandonment thus influences soil nitrogen content and C/N ratio. Near complete recovery of carbon is predicted to take several centuries [92].

The documented relationship between soil carbon and agricultural history, along with present study results indicating a strong correlation between soil carbon and nitrogen (R = 0.94, p < 0.01), suggest an additional confounding factor: protected areas with higher housing density in the surrounding landscape may have gone through an earlier process of agricultural abandonment. Supporting this idea, Wisconsin and Massachusetts data each demonstrate a high correlation between population density and conversion from cropland to grassland [93,94]. While we found no significant relationship between housing density and date of the oldest housing development surrounding each site, this does not preclude a potential correlation between housing density and the more poorly documented date of agricultural abandonment.

Date of agricultural abandonment would also presumably affect the quantity of overstory biomass. The relationship between overstory biomass and date of agricultural abandonment is intuitively due to the positive relationship between total stand biomass and stand age [95], but biomass on old fields, or abandoned agricultural land, is also correlated to soil N-pool size [96]. The possibility that older protected areas in our study region had larger trees and greater overstory biomass would be consistent with the relationship between site age and tree biomass observed in the urbanizing landscape of the Colorado Front Range [17,18]. Our study focused on spatial relationships in the present-day landscapes in Southeastern Michigan, but future additional research could more thoroughly investigate the temporal components of relationships between age of protected areas and the time that these have been in the proximity of housing developments, related to time since conversion from agriculture. Further research could be conducted to determine whether age of protected sites controls overstory biomass and also whether age of protected sites correlates with temporal trends in housing construction, including types of housing developments in the region, their soil and vegetation characteristics, and whether these have changed over time. Possibly, as more of the landscape has been converted to residential land use each decade, a trend toward increasing lot size in some exurban developments could mean that older patches of protected forest are more likely to be surrounded by greater housing densities. Further research would be needed to test this idea.

4. Conclusions

Our results supported the hypothesis that soil and vegetation would be more N-rich in edges of protected areas surrounded by greater housing density. Our results also indicated that in the fragmented urban to rural landscape of Southeastern Michigan, in a weak but marginally significant manner, soil N concentrations in the forest edge zones of conservation reserves were highest where surrounding housing densities were greatest. These findings are similar to weak correlations observed between urban densities and measures of ecosystem N cycling found in other landscape-scale studies in other regions. We found stronger control over soil N concentration associated with landscape position and site physiography, which confound efforts to directly test the effects of surrounding housing density on the biogeochemistry of protected areas.

Nitrogen concentration in foliage, an additional potential indicator of site N availability, did not exhibit strong differences among sites, but did demonstrate an overall moderately rich level compared against forest trees in other regions.

The current study demonstrates that an indirect gradient sensu McDonnell et al. [21,22], with refinement, has the potential to be an appropriate approach for the study of biogeochemistry along an urban–rural continuum. Housing density at the census block group scale was found to be relevant to ecosystem nutrient cycling, both in terms of soil nitrogen content and overstory biomass. Urban effects were thus found to be correlated to the local surroundings, rather than simply correlated to distance from an urban core as in a linear belt transect. The non-linear (loge) relationship between LA/HU and biogeochemical properties, while weak, suggests that the first few housing units in an area are the most impactful. Future studies should include more detailed soil measurements and should refine landscape-scale approaches to address the confounding factors of landscape position and site age.

Acknowledgments

The authors thank field assistant Caitlin Shrestha and the managers of the 12 field sites. Dr. Donald Zak, Rima Upchurch, and others in the SNRE Soil and Plant Laboratory provided valuable assistance. We also thank Stephanie Hart and two anonymous reviewers whose insights improved the manuscript. This research was partially funded by the Rackham Graduate School, the School of Natural Resources and Environment at the University of Michigan, and an award from the National Science Foundation, Dynamics of Coupled Natural and Human Systems Program (award GEO 0814542).

Author Contributions

Ari Kahan, William Currie, and Dan Brown conceived and designed the study. Ari Kahan carried out field measurements. Ari Kahan and William Currie analyzed data. Ari Kahan, William Currie, and Dan Brown prepared the manuscript for publication.