Introduction

Shrubs are encroaching into grasslands globally resulting in altered landscapes and ecosystem functions (Tape et al. 2006; Archer et al. 2017; Ratajczak et al. 2017; Sabater et al. 2020; Scharnagl et al. 2019; Huang et al. 2020). Site-specific drivers include warming temperatures, changes in frequency and amount of precipitation, land-use changes, as well as altered grazing and fire regimes (Huang et al. 2020; Malfasi and Cannone 2020; Sulwiński et al. 2020; Tormo et al. 2020). Climate forecasts predict more variable precipitation patterns, including amount and frequency of rainfall and storm events (Liu et al. 2017). Despite this complexity of influential abiotic factors, shrub expansion into grasslands typically leads to increased annual net primary production (ANPP) and altered ecosystem dynamics (Knapp et al. 2008; Gherardi and Sala 2015; Havrilchak and Schedlbauer 2018; Higginbottom and Symeonakis 2020). For example, water use can be dramatically altered by a shift in dominant plant functional type as demonstrated in arid systems (Medeiros and Pockman 2011; Iqbal et al. 2021). Trees and shrubs generally have higher primary production and thus use more water, thereby altering water availability for surrounding species (Winkler et al. 2019; Peguero-Pina et al. 2020).

Across coastal landscapes, encroaching woody vegetation may act as ecosystem engineers, modifying microclimate to reduce mortality, increase ANPP, and accelerate woody expansion into grassland (D’Odorico et al. 2010; He et al. 2015; Wood et al. 2020), but negatively influence preexisting grassland species. Vegetation growing in sandy soils of barrier islands is particularly vulnerable to changes in water availability due to the reliance on a precipitation-fed freshwater lens in the soil, which temporally varies due to both inundation and drought (Aguilar et al. 2012; Kiflai et al. 2020). Increased nitrogen, phosphorus, and potassium input from leaf litter can accumulate over time and create nutrient rich areas of soil or “islands of fertility” for older woody vegetation (Brantley and Young 2010; McCulley and Jackson 2012; Farella et al. 2020; Zhou et al. 2021). These small-scale changes eventually lead to cumulative landscape change, including increased carbon assimilation and primary production, as well as altered landscape hydraulics and function (Eldridge et al. 2011; Nadal-Romero et al. 2018; Ward et al. 2018; Zinnert et al. 2019). When shifts in climate are coupled with oceanic drivers (e.g., sea-level rise and storm surge/inundation), coastal systems are extremely sensitive to the pressures and pulses of these drivers (Ezer et al. 2013; Pivovaroff et al. 2015; Charles et al. 2020). Barrier islands make up a small portion of global coastlines, but they protect significant areas of economic and ecological value (Stutz and Pilkey 2001; Feagin et al. 2015; Hanley et al. 2020).

Along the US Atlantic and Gulf coasts, woody cover, primarily the native, evergreen shrub, Morella cerifera, has been expanding into coastal grasslands (Battaglia et al. 2007; Zinnert et al. 2019). On the Virginia barrier islands, expansion of this salt sensitive species has occurred despite a decrease in island surface area (Zinnert et al. 2016). The main driver of M. cerifera expansion is warming winter temperatures in the region, 0.06 °C year −1 since the 1950s (Huang et al. 2020; Goldstein et al. 2018). Perhaps contributing to success in sandy coastal soils, M. cerifera has an association with a nitrogen-fixing actinomycete, Frankia, in root nodules (Young 1992; Brantley and Young 2010) that ensures an adequate supply of nitrogen (N), one of the most limiting nutrients in terrestrial coastal ecosystems (Brantley and Young 2010).

Morella establishes in low-lying elevations (i.e., swales) behind foredunes that are ~ 1.5 – 2 m in height, where the shrubs have access to the freshwater lens and are protected from salt spray and intrusion from storm surge (Young et al. 2011; Woods et al. 2019). Adequate access to N via N-fixing actinomycete root nodules allows for high ANPP, in an otherwise nitrogen limited habitat, and development of dense, closed canopy thickets (Shiflett and Young 2010; Shiflett et al. 2013, 2014; Zinnert et al. 2016). Although M. cerifera influences on microclimate and species composition have been quantified on a small scale (Thompson et al. 2017; Wood et al. 2020), effects of annual macroclimatic patterns on primary productivity over longer timeframes have not been documented.

Ecological studies often rely on chronosequences (i.e., space for time substitutions) to investigate succession of plant communities, diversity, primary production, soil nutrients, soil microbial activity, and fungal communities over long timeframes (Springsteen et al. 2010; Hollingsworth et al. 2010; Blaser et al. 2014; Cava et al. 2018; Yang et al. 2019; Li et al. 2020). Chronosequences enable long-term studies in a more compressed time scale. They may occur in areas with regular disturbance or landscape change, such as fire, overwash, flooding, restoration, or sand accretion (Brantley and Young 2010; Hollingsworth et al. 2010; Pellis et al. 2019; Laskar et al. 2021). On Hog Island, Virginia, sand accretion occurred on the north end of the island at ~ 5 m year −1 since the late 1800s, establishing a chronosequence of dune-swale complexes of progressively younger soils from the bayside to the ocean side of the island (Fig. 1) (Hayden et al. 1991; Shiflett et al. 2014). This phenomenon allows us to study differences among shrub and soil ages and drivers of ANPP over time.

Fig. 1
figure 1

Location of Hog Island along the east coast of Virginia. Star denotes the NOAA weather station located in Painter, Virginia. Areal image of the northern end of the Hog Island, Virginia chronosequence with thickets decreasing in age from left to right: bayside shrub (BS), mid-island shrub (MS), young shrub (YS), grassland (GR), colonizing shrub (CS). Location of wells are noted for the BS, MS, and YS shrub thickets with yellow dots

Our study determined the effects of 1) macroclimate changes on Morella cerifera ANPP using long-term data (1990 – 2007) across the Hog Island, Virginia chronosequence and effects of 2) shrub encroachment on soil N and carbon (C), as well as groundwater dissolved organic carbon (DOC) and nitrogen (DON). We hypothesized that 1) temperature and hydrological variables (i.e., water table depth and precipitation) would be influential controls on shrub ANPP due to photosynthetic temperature optimum for M. cerifera and known reliance on the soil freshwater lens (Young et al. 2007; Aguilar et al. 2012). We further hypothesized that 2) soil N and C would be elevated in shrubland compared to grassland and highest under older thickets, due to the greater input of leaf litter through time. We expected groundwater DOC and DON would be more variable due to potential hydrological movement.

Methods

Study Site

Our study was performed on Hog Island, Virginia, part of the network of barrier islands and marshes that comprise the Virginia Coast Reserve (VCR), owned by the Nature Conservancy and supported by the National Science Foundation Long-term Ecological Research (LTER) program. The island is located ~ 10 km east of the Eastern Shore of Virginia, USA (37°270 N, 75°400 W) (Fig. 1). Remote sensing analyses have documented landscape changes, including sand accretion on the north end of Hog Island, as well as subsequent establishment of Morella cerifera shrub thickets (Young et al. 2007; Zinnert et al. 2016). The northern end accreted ~ 5 m yr−1 during the timeframe of this study, resulting in a chronosequence moving from the bayside marsh to the ocean, where soil age increased with distance from the ocean (Young et al. 1995). Through the accretion process, a series of dune ridges formed that protect low-lying elevation swales from salinity and other disturbance (Young et al. 1995, 2011). Over time, the swales converted from grassland to dense shrub thickets dominated by M. cerifera of different ages that varied in stem size-class structure (Brantley and Young 2007; Zinnert et al. 2016). The four thicket ages used in this study from ocean side of the island to mainland side were designated as colonizing shrub (CS, ~ 10 yr), young shrub (YS, ~ 20 yr), mid-island (MS, ~ 40 yr), and bayside (BS, ~ 60 yr). Each thicket was > 10 m wide and extended for > 50 m parallel to the ocean shoreline (Fig. 1). All shrub thickets were established behind dunes which reduce stress from salt spray and blowing sand (Young et al. 2007, 2011).

Primary Production

Beginning in 1990, in mid-October (after growth had ceased but before significant leaf fall) annual aboveground biomass was harvested from the current-year growth, which is differentiated by stem color (Young 2007). Along the chronosequence, fifty shoots were collected in total from five random locations representing different shrubs within each of the four chronosequence thickets. Prior to sample collections shoots were examined to ensure that all originated from different mainstems. Each shoot was dried at 80 ºC for 48 h and weighed. New shoot growth (NSG) was expressed as mass and then input into allometric equations to estimate ANPP as derived from Brantley (2005). Incremental aboveground biomass, i.e., mainstem growth rings were not included. The allometric equations included annual leaf production (ALP) determined from seasonal leaf litter collections and twig mass: leaf mass ratios calculated from separating twigs and leaves from NSG. Equations were produced for each thicket (Brantley and Young 2008). Due to constraints of site access, ANPP for the M. cerifera thickets was collected in all years except 1991.

The general equation for estimating ANPP is:

$$ANPP = {\text{ALP}} (g {m}^{2} {yr}^{-1}) + ((\mathrm{Twig Mass}:\mathrm{LeafMass }* NSG (g))$$

Allometric equations for each thicket are:

$$\begin{array}{c}BS ANPP=310 +(57*NSG)\\ MS ANPP=344+(81*NSG)\\ \begin{array}{c}YS ANPP=360+(108*NSG)\\ CS ANPP=707+(31*NSG)\end{array}\end{array}$$

Climate

We accessed data (1990 – 2007) from the VCR database to identify important climatic variables related to shrub ANPP. Daily precipitation, hourly temperature, and windspeed data were acquired for both Hog Island and Oyster, VA meteorological stations based on data availability (Porter 2023). Water table depth data were obtained from long-term monitoring wells in the BS, MI, and YS shrub thickets on Hog Island, VA from 1990–2007 (Day 2007). Precipitation was quantified as an annual sum and daily mean. Monthly mean sea-level (MSL) data were accessed from NOAA (NOAA 2023). MSL data were summarized to yearly means. Temperature, wind speed, and water table depth were summarized to mean, minima, and maxima for each year. Temperature data were further summarized to the number of days above 29.5 °C to account for the photosynthetic optimum of M. cerifera (Young 1992).

Soil Characteristics

Soils were collected haphazardly from CS, YS, MS, and BS shrub thickets across the chronosequence and also from adjacent grasslands (GR) in 2017 (n = 3) using an aluminum tube inserted to a depth of 50 cm. Soils were separated into 10 cm sections and each section of known volume was dried at 105 °C for 72 h and weighed to determine bulk density. Soil organic matter content was quantified by loss-on-ignition following removal of any root material. Samples (1 g) were baked in a muffle furnace at 550 °C for 5 h and reweighed to estimate soil organic matter content (%). Subsamples from each 10 cm depth were sent to Cornell Isotope Lab (COIL, https://cobsil.cornell.edu/) for analysis of soil C and N.

Groundwater

Groundwater was collected in 2018 in BS, YS, and GR plots (n = 6). Temporary wells (1 m depth) were dug using a soil auger and fitted with PVC pipe. Wells breached the freshwater lens and were allowed to fill with groundwater, which was hand pumped into ~ 275 mL bottles at ~ 50 cm depth. Bottles were kept on ice in the field and frozen until further processing. Groundwater samples were shipped frozen on dry ice and processed by the University of Georgia Agricultural and Environmental Services Lab (http://aesl.ces.uga.edu/) for dissolved organic carbon (DOC) and dissolved organic nitrogen (DON).

Statistics

ANPP, bulk density, and organic matter content were log-transformed due to non-normality. Linear regressions were performed to determine trends in temperature data over time. Analysis of covariance (ANCOVA) was used to assess trends in water table metrics over time and to determine ANPP differences among thickets. No significant interactions between time and thicket were found, so models were run with time as a covariate. Tukey post hoc tests were used to quantify differences among thickets. Pearson correlations determined the relationships between precipitation, MSL, and water table depth by thicket. The influence of environmental drivers over time on ANPP (not transformed) was quantified with generalized linear models (GLMs). GLMs with a lognormal error distribution were used to assess the effects of the following predictor variables on ANPP: thicket age, days above 29.5 °C, annual maximum wind speed, annual mean, maximum and minimum water table depth, mean, maximum and minimum temperature, and annual sum and daily mean precipitation. Best-fit models were identified by AIC corrected for small sample size (AICc, Burnham and Anderson 2002). Due to lack of well data, CS thickets were excluded from GLM analysis. Soil variables were compared among thicket age and soil depth with two-way ANOVA. Significant (P ≤ 0.05) interactions were followed by Tukey post hoc tests. When interactions were not significant, main effects comparisons were analyzed using Tukey post hoc tests. Significant variations in DON and DOC were identified using one-way ANOVAs.

Results

Many of the climatic variables showed high variability but there were clear trends across the timeframe of the study, 1990–2007. Mean air temperature increased through time (β = 0.16, r2 = 0.21, P < 0.001) (Fig. 2). Maximum air temperature increased 0.32 °C annually but minimum temperature did not have a significant trend (P = 0.02 and P = 0.27, respectively). Similarly, mean and minimum water table depth decreased from the soil surface significantly over time (mean: β = 1.13, r2 = 0.53, P < 0.001; minimum: β = 2.92, r2 = 0.34, P < 0.001; maximum: P = 0.12) (Fig. 2). In contrast, annual precipitation was variable over time but did not exhibit a trend (P = 0.14) (Fig. 2). There were no correlations among annual or daily mean precipitation and water table depth metrics (P > 0.05), indicating other influences on water table depth. MSL was correlated with mean water table depth at the BS (r = 0.74, P = 0.0008) and MS locations (r = 0.59, P = 0.012), but not at the YS thicket. Mean and maximum water table depths were furthest from the soil surface in YS (mean: 3.2 cm, maximum: 34.3 cm) and above or closest to the soil surface in BS (mean = -13.7 cm, maximum: 12.2 cm) (mean: F = 17. 01, P < 0.0001, maximum: F = P = 0.014). No differences were observed in minimum water table depth.

Fig. 2
figure 2

Relationship between climatic variables and time from 1990–2007. Annual precipitation (ppt) did not exhibit a signficant relationship over the study timeframe (A), but mean temperature (B), mean water table depth (C), and minumum water table depth (D) had positive signficant relationships over 17 years (α = 0.05)

For the entire study period, younger M. cerifera thickets (CS and YS) had significantly higher ANPP than the MS and BS thickets, with BS exhibiting significantly lower values than MS (Fig. 3). Year and thicket age influenced ANPP (P < 0.001 for both) with no interaction, indicating that as thickets age, production decreased – YS had 1176 g m−2 y−1 higher ANPP than BS in 1990 and by 2007 had 270 g m−2 y−1 higher ANPP than BS. In 2007, CS, the youngest thicket had only 429 g m−2 y-1 higher ANPP than BS (Fig. 3).

Fig. 3
figure 3

Mean ANPP (g m−2 y−1) for all shrub thickets across whole time scale (top). Bayside shrub = BS, mid-island shrub = MS, young shrub = YS, colonizing shrub = CS. Differing letters indicate statistical differences via a post-hoc Tukey test (α = 0.05). ANPP among thickets over time (1990—2007, bottom)

The variation observed in ANPP across all sites over time indicates possible environmental drivers influencing ANPP rather than shrub age being most important, thus we tested several parameters (i.e., air temperature, precipitation, water table depth). The top three models in the accounted for 74–79% of variation and included strong influences by thicket age and hydrological characteristics. The best fit model with the lowest AIC had five variables – thicket (BS-MI), thicket (YS-MI), max water table depth, min water table depth, and mean precipitation (Table 1). ANPP declined with thicket age and higher water table maximum and minimum, whereas mean precipitation increased ANPP. The second best-fit model included maximum temperature but otherwise maintained the same terms as the top model (Table 1).

Table 1 Parameter estimates, Chi-square test statistic, P-value, AICc, and generalized R2 for predictor variables in top GLMs assessing shrub ANPP

Bulk density did not differ among thicket age and grassland or with soil depth (Table 2); however, several soil characteristics were influenced by shrub presence, thicket age, and depth. Organic matter content (%) was higher in older shrub thicket soils (F = 4.8, P = 0.002; Table 2), but did not differ by depth. The oldest soils associated with BS thickets showed increased nutrient accumulation. Soil N content differed by depth (F = 6.2, P < 0.001) and thicket (F = 2.7, P = 0.042), but there was no interaction (Fig. 4). The lowest N content occurred at the 20 – 30 cm and 30 – 40 cm depths; Soil N in GR plots was lowest by 2.7 g m−2 (Fig. 4). Soil C content followed a similar trend to soil N (depth F = 5.9 and P < 0.001, thicket F = 6.4 and P < 0.001); soil C at 20 – 30 was lower than 0 – 10 cm by ~ 260 g m−2 and 30–40 cm was lower by 320 g m−2. C in GR plots was 150 – 216 g m−2 lower than the BS, MS, and YS thickets while CS thickets had intermediate C levels (Fig. 4). There was a gradient in mean values of C and N of oldest thickets to grassland, indicating the time is needed for the high litter input to accumulate, break down, and influence soil values. DON in the water table was higher under BS than YS or GR (Fig. 5). There was high variation in DOC under older thickets, which only differed from DOC in GR but was not statistically different than DOC in YS water table (Fig. 5).

Table 2 Bulk density (g m3) and soil organic matter content (%) among shrub thickets and grassland by depth class
Fig. 4
figure 4

Effect of depth and plot location (BS = bayside shrub, MS = mid-island shrub, YS = young shrub, CS = colonizing shrub, GR = grassland) on soil N (A, B) and soil C (C, D) collected in 2018 (n = 3). Different letters indicate statistically significant differences using a Tukey post-hoc analysis (α = 0.05)

Fig. 5
figure 5

DON and DOC from water table in bayside (BS), young shrub (YS), and grassland (GR) plots collected in 2018 (n = 6). Different letters indicate statistical differences using a post-hoc Tukey test (α = 0.05)

Discussion

Shrub encroachment in coastal systems has local consequences – microclimate modification and changes in soil moisture and community composition (Thompson et al. 2017; Wood et al. 2020; Wood 2021). Our study sought to better understand the landscape-scale impact of shrub expansion on carbon dynamics and to identify meteorological inputs that influence shrub primary production. Shrub encroachment results from a variety of drivers and controls on ANPP, depending on the system and the species (D’Odorico et al. 2012; Huang et al. 2020; Malfasi and Cannone 2020; Sulwiński et al. 2020; Tormo et al. 2020). The chronosequence in our system offers a unique opportunity to better understand potential drivers of shrub ANPP and ecological consequences of encroachment in a mesic system. With warming winters driving shrub cover expansion (Huang et al. 2020) and predicted changes in rainfall, we expected vegetation that is reliant on the precipitation-fed freshwater lens to be sensitive to macroclimatic water dynamics and temperature extremes.

Thicket age exhibited a strong effect on shrub ANPP, supporting previous research that observed decreased physiological activity and leaf area index as shrubs age (Brantley and Young 2007; Shiflett et al. 2014). This contrasts with other systems and species that do not exhibit decreased ANPP with age (Rundel and Parsons 1979; Hughes et al. 2006). Hydrology, both precipitation and depth to the water table, were the other main factors influencing ANPP; higher precipitation resulted in greater ANPP, whereas water table depth became closer to the surface over time and was associated with a small decrease in ANPP. Mesic soils typically favor shrub establishment over grasses (Gremer et al. 2018) and in our system, shrubs establish only in low-lying swales where water availability is closer to the soil surface (Young 1992), but they are sensitive to flooding (Tolliver et al. 1997). Precipitation and water table controls on ANPP have been observed via normalized difference vegetation index (NDVI) in our system (Aguilar et al. 2012). During wet years, shrub NDVI (a surrogate for primary production) was correlated with precipitation, whereas in dry years, NDVI is correlated with depth to the water table. These relationships depend on shrub age and in the coastal environment, are controlled by rising sea-level.

Sea-level rise is high in the mid-Atlantic (5.9 mm y−1), with evidence that it is accelerating (Sallenger et al. 2012). At our site, yearly mean sea-level (MSL) was correlated with mean water depth in the two older shrub thickets. Lack of correlation with the water table under the youngest, most productive shrub thicket was likely due to the high leaf area index, hydraulic conductivity, and transpiration by M. cerifera (Brantley and Young 2007; Shiflett et al. 2017; Wood et al. 2020). With projected variability in precipitation altering water table recharge and potential salt intrusion by more frequent and intense storms, the water table will likely continue to play an important role in landscape ANPP (Groeneveld 2008; Bryan et al. 2016; Geißler et al. 2019; Kiflai et al. 2020; Holdrege et al. 2021). As sea-level rises and precipitation amount and frequency shifts with climate change, the impact on ANPP of woody vegetation in mesic environments will likely intensify (Liu et al. 2017; Xu et al. 2018).

Temperature plays a role in expansion and establishment of shrub cover in many ecosystems, and macroclimatic warming since 1950s is attributed to M. cerifera expansion in our system (D’Odorico et al. 2010; Osland et al. 2013; Gabler et al. 2017; Huang et al. 2020), but it did not correlate with ANPP over a 17-year timescale. Morella cerifera primary production is less sensitive to extreme temperatures perhaps because of microclimate modification (Thompson et al. 2017; Sinclair et al. 2020; Wood et al. 2020). The decrease in extreme temperatures by M. cerifera structure and physiology keeps the temperature closer to the photosynthetic optimum of 30 ºC (Young 1992). Our results contrast with other systems where increased shrub primary production is linked to warming temperature and less impacted by water availability (Filella et al. 2004).

Shrub establishment altered the soil characteristics are altered over time; the oldest soils (BS and MI) had higher N and C content than younger soils (CS and GR) indicating that N and C accumulates over time (Raciti et al. 2011; Lunstrum and Chen 2014; Throop and Lajtha 2018; Lan et al. 2021). Similar findings were observed in soil organic matter content, although it did not increase with soil depth. Despite this accumulation, it is still markedly lower when compared to forested habitats (Bell et al. 2021). However, DON and DOC were much higher in the oldest thicket, indicating leaching into groundwater, thus, soil N, C, and organic matter content may not be reliable indicators in sandy systems. As N and C accumulate in the freshwater lens, there is an increased likelihood that the nutrients may move throughout the water table and enrich other parts of the landscape (Hefting et al. 2004; Muwamba et al. 2018). This may be especially important in the relatively nutrient poor, sandy soils on barrier islands. Some mesic habitats vary in woody expansion influence on N and C; soil C is lost due to increased microbial activity and N is accumulated from litter input (McCulley and Jackson 2012). Other mesic grasslands exhibit notable increase in belowground C and no change in N with older woody vegetation (Ward et al. 2018; Charles et al. 2020), partially due to decreased soil biomass breakdown rates. There was a depression in both soil C and N at the 20 cm – 40 cm depth. This soil zone could experience variations in water table saturation that could facilitate more nutrient leaching and an altered microbial community (Masterson et al. 2014; Hashim et al. 2019; Xia et al. 2021); more research is needed to understand declines at these depths.

Although mesic shrub establishment in our system is largely driven by warming temperatures (Huang et al. 2020), ANPP was principally associated with shrub age and water availability. The rate at which shrubs decrease ANPP with age is similar regardless of shrub age or year. Water availability influence on mesic ANPP is understudied but plays a significant role; precipitation increases ANPP while higher water table decreases ANPP. Water table depth is influenced by SLR but may also modified by water use in younger shrubs (Shiflett et al. 2017), further complicating the relationship between shrub presence and water table depth. The soil nutrient dynamics of shrub encroachment are complex because of nutrient leaching and the variability of water table depth over time. In general, soil N and C accumulate with shrub age, but leaching into groundwater occurs and the fate in groundwater is unknown (Masterson et al. 2014). Coastal habitats have additional soil nutrient flux due to loss of N and C from the habitat via erosion (Martínez et al. 2019); sediment movement through erosion is increased by shrub establishment in our system (Reeves et al. 2022). As shrubs are globally expanding into suitable habitat at unprecedented rates (Zinnert et al. 2019; Huang et al. 2020) these environmental alterations have big impacts in in carbon assimilation and water dynamics of mesic systems.

Conclusion

Due to the paucity of studies in mesic coastal environments, our study sought to better understand the landscape-scale impact of shrub expansion on carbon dynamics and to identify meteorological inputs that influence shrub ANPP across a barrier island landscape. Our results contrast with other systems where increased shrub ANPP is linked to warming temperature and less impacted by water availability. Environmental drivers influencing the encroachment of shrubs in this coastal system (i.e., warming temperature) are not important in predicting shrub ANPP. Rather, complex interactions between shrub age and hydrological properties impact primary production, contributing to carbon storage in barrier island soils. Given the dynamic nature of coastal systems and sediment movement, landscape predictions of coastal carbon should include age of soils, vegetation, and hydrology of the system.