Work Description

Date: 22 August, 2021

Dataset Title: Collapse of a Black Tern Colony (Chlidonias niger) as a Result of Climate Change Driven Lake-Level Extremes and Anthropogenic Habitat Alteration

Dataset Creators: J.L. Fuller, E. Rowan, A. Landgraf, K. Alofs, J. Foufopoulos, A. Gronewold

Dataset Contact: Jennifer Fuller jenful@umich.edu

Funding: (SEAS Master's Thesis Grant), (Rackham Graduate Student Research Grant), (The E. Alexander Bergstrom Research Award), R/CGLH-12 (The Michigan Sea Grant Graduate Student Fellowship, NOAA)

Key Points:
- We analysed 8 years (2013-2020) of black tern breeding success data, lake levels, and remote-sensing derived habitat data specific to each nest
- We applied a GLM (General Linear Model) to capture nest success regarding incubation timing and habitat variables. We applied the result of the GLM's habitat relationships to ArcMap layers to visualize survival likelihood on a spatial scale. We also generated a GLMM (General Linear Mixed Model) to calculate the relationship between sub-colony population size and the amount of available nesting habitat.
- Black tern breeding success significantly coincided with shallower water levels, further distance from the main lake, and less surrounding dense vegetation.
- Black tern breeding success and population size both dropped considerably between 2013 and 2020, coinciding with record lake level rise and apparent depletion of safe breeding habitat

Research Overview:
Global climate change is expected to interact with existing environmental stressors and increasingly impact biodiversity. Great Lakes coastal wetlands and wildlife are under multiple threats including development, invasive species, and pollutants. How such stressors interact with climate-change initiated, extreme lake level fluctuations and alter wildlife populations is not well understood. We investigate the rapid decline of black tern (Chlidonias niger) over eight years in one of Michigan’s largest colonies, Lake St. Clair. Using field monitoring and remotely sensed data, we found ideal breeding habitat could not shift upland, with rising water, as coastlines were developed or invaded by Phragmites australis. Subjected to progressively deeper and unstable habitat, nests were likely more exposed to inclement weather and aquatic predators. We present a compelling case demonstrating how climate change and multiple stressors are major conservation management concerns for black terns and wetland biodiversity.

Methodology:

1. Study Region

As the largest inland freshwater delta in the world, St. Clair Flats covers roughly 101 km2 stretched into the 1,114 km2 Lake St. Clair. St. Clair Flats is likely a major attraction for black terns given their preference for expansive wetlands alongside bodies of open water (31, 45, 54). St. Clair Flats is historically home to one of the largest colonies of black terns in the Great Lakes Region, consisting of 137-400 breeding pairs. The capacity of St. Clair Flats to host such a large population of black terns and other species gives the region high conservation value, also reflected in its designation as a globally Important Bird Area (IBA) (17).
The coastal freshwater wetlands are generally dominated by native Schoenoplectus sp. (bulrush), Typha spp. (cattail), and an invasive clonal dominant, Phragmites australis (common reed). Invasions have been facilitated by a combination of heightened anthropogenic pressures; furthermore, invasive species establishment was likely facilitated by a period of low lake-levels (1999-2013) (48, 49). Black terns build nests almost exclusively from the broken stems of emergent vegetation. Most nests sit on top of floating aggregations, or mats, of these stems, though they also frequently take advantage of logs, wood planks, or floating pieces of Styrofoam to use as a platform. In other regions, black terns also use shallow sedge tussock habitat, but there do not appear to be areas of this type of habitat supporting the St. Clair population.

2. Hatch Success Materials

Tern breeding colonies at St. Clair Flats (SCF) (2013-2020) were monitored by no less than two volunteer and staff research technicians at least 1-2 times a week between spring migration/arrival (~May 15) and fall migration (~July 30). Monitoring colonies that were small and isolated were at times given lower priority and visited less frequently. Population sizes were estimated using records of head counts taken throughout the season at each sub-colony. Details of this procedure are in Supplementary Materials (S1).
Sub-colonies were initially flushed to estimate the number of nesting pairs in the area. Nests were then located by pinpointing where the adults landed after flushing and were subsequently georeferenced using a handheld GPS (2013-2016) or ArcGIS Collector App (2017-current). Researchers collected data on the dominant vegetation type(s) used for nest-construction, and on water depth measurements at each nest (when possible) using a marked PVC pole (2018-2020). The number of nests within a 30-m radius were quantified post-field season using the “near” function in ArcMap amongst nest GPS points from the same year.
Given black tern re-nesting ability, colonial nature, and high chick mobility, identifying nest survival required careful observations by field researchers and an understanding of the nest’s age when possible. Except for cold-condition visits early in the season, eggs were given a “float” test to estimate clutch initiation (egg age), hatch date, and prioritize revisits for capture and banding. Records of these tests or estimated hatch dates were less consistent in early years, but this improved over time. Weights from chick banding data was secondarily used to estimate their age followed by the age of the nest based on an average incubation time frame of 21-22 days. The method with the greatest accuracy was chosen for estimating a nest age if both a float test and chick weight occurred (e.g., the nest was found the day the clutch initiated, and therefore more accurate than using the chick’s age-weight). Additional information on how clutch initiation was calculated is found in the Supplementary Materials (S2).
Each nest was visited two to four times during the field season to determine the nest’s outcome. For this work, categories of nest outcomes were combined to determine whether nest eggs either ‘failed’ (n = 165) or ‘hatched’ (n = 286, regardless of chick survival and fledging), leaving 370 nests with unknown final status See Supplementary Materials (S3) for a breakdown of nest outcomes. The 95% confidence interval for each year’s hatching success rate was calculated with a beta distribution using the qbeta() function from the ‘ExtDist’ package in R (55).

3. Water Depth Materials

Water depth was collected and calculated using a combination of high spatial resolution bathymetry (3 m) from the National Oceanic and Atmospheric Administration (NOAA) Office for Coastal Management Lake Level Viewer Dataset and surface water levels from the NOAA Huron Erie Connecting Waterways Forecasting System (HECWFS) model (9). The bathymetric elevation (m) was determined for each nest in ArcMap using the “sample” function from the open source Lake Level Viewer Tool bathymetric map (56). Water depths were then estimated by subtracting the bathymetric elevation from the HECWFS surface elevation, standardized to the season onset (May 15th). Estimated water deaths were standardized to centimeters with zero set as the shallowest relative depth. The methodology for extracting the data required for this calculation is elaborated in the Supplementary Materials (S4, S5).

4. Geospatial Habitat Variable Materials

We chose geospatial habitat variables and sources based on previous literature and the resolutions of available imagery. Black terns are frequently reported to prefer habitats with a balance between open water and vegetation, vegetation percentages ranging between 25-75% (16, 22, 23, 45, 48, 49). Though from a larger spatial scale this measure is useful, it does not fully address the true habitat complexity of their breeding grounds. Not only do vegetation type and structure vary across different wetlands used by black terns, they also can vary significantly within a colony. While monitoring nests in the field in 2019 and 2020, researchers found that larger colonies were usually found in areas that contained plentiful floating, dead plant material, often secured by semi-dispersed Typha spp. or Schoenoplectus sp. Black terns avoided nestin in densely packed Typha spp. and the highly invasive Phragmites australis which impede takeoff and landing, visibility, and prevent mat build-up required for nest construction. We therefore used remote sensing methods to reconstruct vegetation classes corresponding to nesting data collected over previous years. Ground truthing was used to determine that habitat structure could be generalized into four classes based on stand density: 1) dense, standing vegetation, 2) mat and scattered or cut vegetation, 3) sparse vegetation and/or sparse, scattered mat, 4) open water. Another important consideration in terms of nesting habitat is scale. The scale of habitat structure impacting the nest could be very fine, as nests are generally less than 12 inches in diameter (58). Previous studies have examined vegetation within a 12-m radius (57) or 2-m radius of nests (31), lower than available remote imagery spatial resolutions. We compared GLM p-values of hatch success in response to 3-m, a median value at 7-m, and 12-m and found that 7-m generated the strongest model results, therefore it was chosen for the final analysis.
To address the need for high-spatial resolution imagery (ideally < 1 m) covering eight years of nest monitoring, we used a combination of open-source aerial photographs from NAIP (National Agriculture Imagery Program) and purchased 4-band commercial satellite images. NDWI (normalized difference water index), which uses green and near infrared wavelengths to delineate water bodies, was chosen to capture the extent of open water (59). NDVI (normalized difference vegetation index) was chosen for capturing average and classified vegetation density estimates as it uses red and near infrared wavelengths to measure photosynthetic concentrations, or “greenness” (60, 61).
We collected 4-band images under 1-m resolution from commercial satellites (Kompsat-2, Triplesat-3, WorldView-2 and 3), and the National Agriculture Imagery Program (NAIP). We chose 1-m resolution images within nine days of each other (standard deviation = 8.71 days) between late June and July to prevent as much timing differences during annual growth seasons as possible. We obtained 5-m resolution imagery (Rapideye-5 and PlanetScope) through Planet Images to capture more general, average NDVI values. With the advantage of higher temporal resolution, images from Planet could be collected for each year during the breeding season, allowing anniversary dates to be interpolated. All images were resampled to 1- or 5-m resolution depending on their source, and geometrically, radiometrically, and atmospherically corrected based on their individual requirements. Imagery preprocessing methods, image dates, and resolutions are detailed in the Supplementary Materials (S6).
Mean NDVI was sampled from each nest’s 7-m radius buffer region from linearly interpolated anniversary images (PlanetScope and Rapideye-5, 5-m resolution). Extraction used the Zonal Statistics 2 toolbox which is capable of handling overlapping polygons. Open water was classified using Natural Breaks (Jenks) Unsupervised Classification on NDWI raster images from yearly 1-m resolution images. Island developments or regions not considered wetland were heads-up digitized using 1-m resolution NAIP imagery from 2014, 2016, and 2018. The open water and island development classes were used to mask 1-m NDVI raster images and generate three vegetation classes using the Iso Cluster Unsupervised Classification tool in ArcMap. The percentage of each class was then sampled within each 7-m radius nest buffer using the Tabulate Intersection Tool. Because there was no available high-resolution imagery in 2015, the percentage of habitat classes surrounding nests were estimated by averaging measured values from classified 2014 and 2016 maps. To measure potential impacts of proximity to the housing developments, any open water, and the larger body of open water (Lake St. Clair), distance values were extracted from the digitized developed regions, the open water class, and a “main lake” delineation using the “near” function in ArcMap. The “main lake” variable was created for each year by delineating the outer edge of the open water class to remove any inundated regions surrounded completely by vegetation. The Supplementary Materials provide further information on the methodology for the “main lake” classification (S7).

5. Hatching Success and Habitat Analysis

The relationships between hatching success and habitat variables were analyzed using R version 4.0.3 (55). A multiple logistic regression with a binomial (logit-link) fit using the glm() function from the ‘stats’ package in R (55) to determine the influence of habitat and biological predictors on nesting survival. The following 10 habitat or biological variables were examined in the analysis as fixed effects: (1) relative, initial water depth (cm), (2) distance to housing developments (m), (3) distance to the main lake (m), (4) distance to any open water (m), (5) percentage of open water within a 7-m radius, (6) percentage of dense vegetation within a 7-m radius, (7) percentage of medium vegetation within a 7-m radius, (8) percentage of sparse vegetation within a 7-m radius, (8) estimated start of incubation as the number of days before or after May 15th of each year, (9) mean NDVI within a 7-m radius, (10) number of nests within a 30-m radius. Continuous habitat variables were first checked for multicollinearity using a Pearson correlation test using the cor() function from the ‘DescTools’ package in R (55). No rs values were greater than 0.6 and all predictors were retained.
To compare model coefficients, all continuous independent variables were normalized using the scale() function of the ‘base’ package in R (55) which computes a z-score for each variable using its mean and standard deviation. Models were compared using stepwise selection in both directions using the step() function from the ‘stats’ package c To determine the best explanation of the data variation, each model was assessed for the lowest Akaike’s information criterion (AIC) (62). A Receiver Operator Curve (ROC) and AUC (area under the curve) was generated using the roc() function from the ‘pROC’ package in R c to measure the best model’s predictive performance (62).

6. Yearly Changes in Lake Level and Habitat Analysis

To visualize the geographic extent of the habitat’s potential to predispose nests’ vulnerability, the coefficients and intercept from our selected GLM were applied to raster layers in ArcMap using the Raster Calculator Tool (see Supplementary Materials S8 for background methodology). For spatial projections, the biological variable describing clutch initiation date was removed. The remaining geospatial variables and associated raster layers were rescaled using min-max feature scaling. The raster calculations generated a final model describing nest failure probability on a scale of 0 to 1. This model’s performance was also assessed using an ROC (Receiver Operator Curve) plot (52) and AUC (area under the curve) using the roc() function from the ‘pROC’ package in R (55).
The final maps were then masked to exclude regions that were determined a priori to be unsuitable for nesting. This includes developed or dry-land islands and peninsulas, open water, and areas where vegetation within a 7-m radius breached an NDVI value of 0.72. Specifically, 0.72 was the maximum NDVI surrounding a known nest and it is assumed that the likelihood of any nest existing within areas with any higher value is extremely low. This was also confirmed by extensive nest searching in the field, which determined that nests are not built within dense monocultures of Phragmites australis or Typha spp. This is because the vegetation mats required for nesting do not accumulate among stands when the vegetation grows too closely together; furthermore, nesting birds are unable to take off or land.
To quantify sub-colony breeding pair population size and their response to changes in sub-colony habitat in the geospatial model, we applied a general linear mixed model (GLMM) using the lmer() function from the ‘lme4’ package in R (55). Predictor variables were chosen a priori, and included the area of open water, uninhabitable vegetation (NDVI>0.72), any habitable area, and area with >50% hatch success. The area of classified predictor variables per sub-colony were extracted from the geospatial model outputs in ArcMap. The response variable, i.e., the number of maximum breeding pairs per sub-colony, as well as the predictor variables were scaled using the scale() function from the ‘base’ package in R (55) prior to running the GLM, to account for considerable differences in measurement units. The area with >50% hatch probability and any habitable area were correlated, as was uninhabitable vegetation and open water extent. The selected model was chosen based on having the highest R2 and including the most significant predictors. After evaluating all possible variable combinations, the selected model included the area with >50% hatch probability and uninhabitable vegetation.

Instrument and/or Software specifications: Remote-sensing habitat data was derived using ArcMap 10.6.1 and data was analyzed in R version 4.0.3

Files contained here:

BLTE_data Folder: R scripts for running full (non-geospatial) analysis, datasets used in R, and figure outputs

1. Code subfolder: R scripts for running full (non-geospatial) analysis

- Word document version of R script. Filename prefix: BLTECode_2021.docx

- PDF version of R script. Filename prefix: BLTECode_2021.pdf

2. data subfolder: Datasets used for nest success GLM analysis, population figure, and nest success x geospatial extent GLMM analysis

- Primary dataset used for the generation of the binomial general linear model to determine habitat and biological variables associated with nesting success in black terns. Filename prefix: nestbeta.csv

- Columns:
1. nestid = Year + Subcolony Code Name + Nest Number
2. subcolony = Code name to designate colony region the nest was located in
3. av.year.wd = Average yearly surface water depth (m) between May and July at Lake St. Clair Flats
4. egg.surv.fl = egg survival code, 0 = failed, 1 = hatched, NA = unknown
5. egg.surv = egg survival alternate code, 1 = failed, 0 = hatched, NA = unknown
6. egg.name = egg survival translated into text
7. start.est = estimated clutch initiation date as the number of days following May 15th of the nests’ corresponding year
8.dem = digital elevation model value (m) overlapping with nest location
9. init.wd = water depth at the onset of the season (May 15th) at each nest location. Values are relative to 0, which represents the nest at the shallowest estimated depth (cm)
10.house.dist = distance between the nest and nearest developed area (m)
11. lake.dist = distance between the nest and the nearest point considered the “main lake” (m)
12. ow.dist = distance between the nest and the nearest point considered any open patch of water (m)
13. thinv7 = percentage of sparce vegetation within a 7m radius of the nest
14. medv7 = percentage of medium vegetation within a 7m radius of the nest
15. thickv7 = percentage of dense vegetation within a 7m radius of the nest
16. water7 = percentage of open water within a 7m radius of the nest
17. ndvi.mean7 = NDVI value (-1 to 1) surrounding a nest within 7m radius
18. ncount30 = number of other black tern nests within a 30m radius of a nest

- Dataset for GLMM, derived from area of different habitat types for used to investigate the relationship between habitat availability and population size. Filename prefix: popgeo.csv

- Columns:
1. subcolony = code name to designate colony region the nest was located in
2. year
3. code = subcolony code + year
4. lake_level = average yearly surface water depth (m) between May and July at Lake St. Clair Flats
5. rawpop = raw population size by subcolony and year
6. area50 = area of habitat within the subcolony with a hatch success probability >50%
7. areahab = area of habitable nesting area
8. areaunhabveg = area of uninhabitable vegetated area
9. areawater = area of open (uninhabitable) water

- Dataset for yearly population changes plotted in Figure 2

- Columns:
1. year
2. population = total estimated number of black terns in the St. Clair Flats region

3. figures subfolder: location for figures 1-6

- Caption: Graphic representation of the different vegetation zones of black tern nesting habitat in order of increasing proximity to land: (A) water lily (Nuphar spp., Nymphea spp.), (B) bulrush (Scirpus spp., Schoenoplectus sp.), (C) cattail (Typha spp.), and (D) sedge (Carex spp.) dominant. As lake levels drop, plant communities migrate lakeward unless propagation is blocked by high hydraulic energy, boat traffic, and/or deep water (E). During drought periods and migration, shallow and exposed soil that no longer supports natural vegetation can become invaded by invasive species Phragmites australis and Typha x glauca (F, red). As lake levels rise again, communities will move up the shoreline, but can be blocked by Phragmites australis establishments or shoreline developments (G). Filename prefix: figure1.pdf

- Caption: Yearly estimated number of black tern breeding pairs across St. Clair Flats (solid red line, right y-axis), average surface lake level (May-July) (m), and yearly hatching success rates (% hatched of total nests and 95% confidence intervals (CI), dots and whiskers, left y-axis). Population size is a “best” estimate of individuals, with maximum and minimums (shown by shaded sections) recorded in 2013 and 2014. Filename prefix: figure2.tif

- Caption: Hatching success probability and 95% confidence intervals with respect to individual variables (A-F) included in the two top binomial general linear models. These include: A. Clutch initiation day (# days after May 15th), B. Percentage of surrounding dense vegetation in a 7-m radius, C. Relative nest water depth (cm), D. Distance to the open water of the main lake (m), E. Distance to developed land (housing) (m), F. and the average NDVI in a 7-m radius. Filename prefix: figure3.pdf

- Caption: Boxplots of the nest variables included in the final GLM, against the respective average May-July surface lake level. Lake levels increased nearly monotonically over the course of the study. The average median for values from 2013-2016 (red) and 2017-2020 (blue) are plotted as lines for reference with respective 95% confidence intervals. Filename prefix: figure4.pdf

- Caption: Decline in tern habitat nesting quality in the study area over the years. Mapped nesting quality is based on the applied binomial GLM output for hatching success as based on geospatial habitat variables (Water depth, Distance to lake & housing, % of surrounding dense vegetation, and mean NDVI). All study years are shown except 2013 (no nest searches) and 2015 (no available imagery). Only those areas for which quality imagery was available for all years have been mapped in detail (inside of two polygons). Water is shown in shades of blue, while land is shown in shades of tan. Emergent wetland, the tern’s nesting habitat, is shown in shades of green ranging from light yellow-green (associated with high nesting success) to green-blue (low nesting success). Nests are shown as red circles. Images show examples of two of the 13 sub-colonies in the study area; Bruckner’s (left) and Canoe (right). Bruckner’s sub-colony was not searched in 2014. While rising lake levels lead to a pronounced reduction of emergent wetlands, shorelines hardened through human activities (roads, dams, build structures) prevent the establishment of new shoreward wetland areas. Filename prefix: figure5.jpg

- Caption: General linear mixed model (GLMM) plots predicting the response of breeding pairs by changes in area of 2 major wetland classes: >50% hatch probability (A) and uninhabitable vegetated regions (B). Colored points denote scaled number of breeding pairs versus the scaled area of each class. Model estimates and 95% confidence intervals are represented by the solid line and filled area, consecutively. Filename prefix: figure6.pdf

BLTE_MasterGIS Folder: All processed raster and polygon datafiles used for the geospatial GLM model, including geospatial output, and nest habitat variables.

1. geomod_inputs subfolder: Raster files used for running GLM model in ArcMap

- Nested subfolders titled by year 2013-2020
- initial water depth on May 15. Filename prefix: initwd(last 2 digits of year, e.g. 13) (raster)
- euclidean lake distance continuous raster (m). Filename prefix: lakedist(last 2 digits of year, e.g. 13) (raster)
- ndvi (input raster, reclassified). Filename prefix: ndvi(last 2 digits of year, e.g. 13) (raster)
- percent of dense vegetation in the surrounding 7-m radius using focal statistics. Filename prefix: thickv(last 2 digits of year, e.g. 13) (raster)
- same euclidean distance from housing used for all years. Filename prefix: housedist (raster)

2. geomod_outputs subfolder: Raster and polygon outputs after geospatial analysis

- comparable clipping extent of each sub-colony for each year. Filename prefix: area_extent (polygon)
- classified_models subfolder
- model values clipped to overlapping regions (1 = lowest chance of failure to 100 = highest chance of failure), with overlaid open water class (150) and uninhabitable dense vegetation (200). Filename prefix: allmosaic(last 2 digits of year, e.g. 13) (raster)
- mixed_models subfolder
- polygon (gridcoded) values for mixed model categories (1 = habitable vegetation, 2 = uninhabitable open water, 3 = uninhabitable vegetation). Filename prefix: alltest(last 2 digits of year, e.g. 13)cl (polygon)
- polygon of regions with only >50% survivability. Filename prefix: only50plus (last 2 digits of year, e.g. 13) (polygon)
- raw_models subfolder
- raw model output after input in arcmap raster calculator tool applying geomod inputs into the GLM. Filename prefix: model(last 2 digits of year, e.g. 13) (raster)

3. nest_habitat subfolder: classifications, NDVI, and distances used for determining nest habitat characteristics

- 1mclasses subfolder
- 1-m classification of open water, thin, medium, and dense vegetation using high resolution NAIP and commercial satellite imagery. Filename prefix: Ks(last 2 digits of year, e.g. 13)all (raster)
- Lakedist subfolder
- 1-m classification of the “main lake” or the St. Clair Flats river delta which does not contain any vegetation. Filename prefix: Lakerasrc(last 2 digits of year, e.g. 13) (raster)
- NDVI_planet subfolder
- 5-m NDVI raster images for each year. Filename prefix: Ndviintmm(last 2 digits of year, e.g. 13) (raster)

Related publication(s):
Fuller, J.F., et al. (2021). Shorebird colony collapses under climate driven lake-level rise and anthropogenic stressors. Forthcoming.

Use and Access:
This data set is made available under Attribution 4.0 International (CC BY 4.0).

To Cite Data:
Fuller, J., Rowan, E., Landgraf, A., Alofs, K., Foufopoulos, J., Gronewold, A. Dataset for "Collapse of a Black Tern Colony (Chlidonias niger) as a Result of Climate Change Driven Lake-Level Extremes and Anthropogenic Habitat Alteration" [Data set], University of Michigan - Deep Blue Data. https://doi.org/10.7302/d104-mg64