This preliminary exploration of marine lichenized fungi (lichens) as bioindicators of water pollution examined the distribution of intertidal lichen communities in the Boston Harbor Islands National Recreation Area with respect to recorded pollution throughout the harbor. We found significant negative associations between pollution measurements and the health of the lichen community based on cover and species richness. We also observed significant differences in species composition between areas of higher pollution and areas of lower pollution, though not enough data are available to establish the pollution sensitivity or tolerance of individual species. We note that difficulties in the collection and identification of marine lichens hamper efforts to use them broadly as bioindicators. This study suggests that marine lichens could prove useful as bioindicators, but more research is needed to understand the differential effects of pollution on individual species as well as to establish practical procedures both for quantifying marine lichen community health and for widespread bioindication using marine lichens. Finally, one species collected during this study, Verrucaria ceuthocarpa, represents a first report for the Boston Harbor Islands National Recreation Area.

We tested potential for the use of marine lichens as bioindicators by examining their distribution in the Boston Harbor Islands NRA and compared those findings with land use histories and water pollution analyses in select locations throughout the park. Six pollution types were examined for this study, data for most of which came from the recent MWRA Boston Harbor Water Quality Updates (2018). The data we used included measurements of nitrogen concentration, chlorophyll-a concentration, total suspended solids, turbidity, qualitative classifications of depositional, erosional, and intermediate sediment types, and the nature of pollution of the Weir River Area of Critical Environmental Concern (Lefebvre et al. 2002; Pahlevan et al. 2018; Taylor 2018).

Collecting procedures were approved by the National Park Service under permit #BOHA-2019-SCI-0005. Six areas were sampled based on accessibility and distance from historical or current point sources of pollution in the Boston Harbor: Bumpkin Island, Grape Island, Lovells Island, Peddocks Island, Thompson Island, and World’s End peninsula. The coastline was also examined at Rowes Wharf in Boston. All sampling occurred throughout the month of August in 2019. We did not take neap and spring tides into account due to difficulties in scheduling travel to the islands; however, the range of the height of the tide on sampling days was approximately two feet (US Harbors 2020). We sampled along the entire length of transects of 1 m width that stretched from the low tide mark to the high tide mark based on mean lower-low and mean higher-high water data (NOAA 2011). In total, 33 transect sites were surveyed. Different islands had different numbers of transects based on the estimated length of suitable rocky habitat along the shore, with each estimated 100 meters corresponding to one sampling site. Satellite images of the Boston Harbor were used to estimate shoreline lengths and to observe possible habitats (Terrametrics 2020). Sites were chosen ahead of time by randomly selecting an integer between 1 and 5000 and estimating the site location by moving that many meters counterclockwise from the landing dock on each island and from the entrance at World’s End. When selected locations did not have suitable habitat for marine lichens, new points were chosen using the same method.

Specimens were identified through microscopy following Taylor (1982). For five specimens that proved difficult to verify morphologically, identifications were confirmed using deoxyribonucleic acid (DNA) extractions and polymerase chain reaction (PCR) amplification of the internal transcribed spacer (ITS) ribosomal DNA region. We used the DNeasy PowerPlant Pro Kit (Qiagen, Valencia, CA) for DNA isolation following the manufacturer’s protocols. Amplification of the ITS was done in 25 μL reactions with 12.5 μL of Promega 2× PCR Master Mix (Promega Co., Madison, WI), 1.25 μL of forward and reverse 10 μM primer, 9.0 μL of H₂O, and 1.0 μL of template DNA. Different primer combinations were tested for optimal results, including ITS1F (Gardes and Bruns 1993) and ITS9mun (Egger 1995) as forward primers, and ITS4 (White et al. 1990) and ITS4A (Larena et al. 1999) as reverse primers. Thermocycler conditions were as follows: initial denaturation at 94 °C for 3:00 min; 35 cycles of denaturation at 94 °C for 1:00 min, annealing at 50 °C for 0:45 min, and extension at 72 °C for 1:30 min; and final extension at 72 °C for 10:00 min (Haelewaters et al. 2018). Purification and Sanger sequencing were outsourced to Genewiz (South Plainfield, NJ). Forward and reverse reads were assembled and edited in Sequencher 5.2.3 software (Gene Codes Co., Ann Arbor, MI), and final sequences were submitted to the National Center for Biotechnology Information  GenBank sequence database (accession numbers MT809480–MT809484). Collections verified with sequences are specified in Table 1. Finally, we compared our samples with voucher specimens from LaGreca et al. (2005). Nomenclature follows Index Fungorum (2020). Specimens are preserved at FH (Farlow Herbarium, Harvard University, Cambridge, MA). Note that field identifications on Bumpkin Island were made using Taylor (1982) and comparison with herbarium specimens from other islands. No voucher specimens were collected from this island  due to practical limitations.

We measured species richness (i.e., the number of species) in each transect and total coverage, size, number of thalli, and reproductive potential (density of perithecia) for each species. Size, cover, and count of lichen thalli within each sampling site were calculated with a combination of measurement in the field (i.e., counts, measurements, and field identifications) and analysis of photographs taken at each site. Size was measured by the longest diameter of each thallus, and cover was estimated by placing a grid on the photograph of each transect and using field notes and pictures to estimate percent cover of the transect. We measured perithecial density with collected specimens by finding the average number of perithecia within three grid templates of 1 cm² placed on each collection using a dissecting

Because our study design and data violated conditions for most parametric statistical analyses (according to the Anderson-Darling test and D’Agostino K² tests for normality), we used nonparametric models and significance tests to analyze our data. For all significance tests, ɑ = 0.05 unless otherwise specified or Bonferroni corrected. We used Spearman’s Rho and associated two-tailed significance tests to assess nonparametric correlation (MacFarland and Yates 2016), with the following hypotheses: H₀: ρ = 0, H_a: ρ ≠ 0. To compare species richness among different areas of differing levels of pollution, we used the Mann-Whitney U test (pairwise) (MacFarland and Yates 2016) and Bonferroni-corrected significance tests with hypotheses H₀: η₁ = η₂, H_a: η₁ ≠ η₂. To show differences in species composition among locations of different levels of pollution, we used non-metric multidimensional scaling (NMDS) using the Bray-Curtis dissimilarity index (McCune and Grace 2002). We used analysis of similarities to test for significant differences between species composition of different groups, again using Bray-Curtis dissimilarity and with the following hypotheses: H₀: there is no difference between the means of two or more groups of (ranked) dissimilarities, and H_a: there is a difference between the means of two or more groups of (ranked) dissimilarities. Finally, to analyze the main contributors to differences in species composition we used the similarity percentage analysis (SIMPER). All tests and models were calculated with the ‘vegan’ package in the R language and environment for statistical computing (Oksanen et al. 2019; R Core Team, Vienna, Austria), with the exception of the pairwise ANOSIM, for which we used PAST software (Oslo, Norway), and SIMPER, for which we used PRIMER-7 (Quest Research Limited, Auckland, New Zealand).

Two README files are attached along with the raw data:

README.lichens_as_marine_bioindicators_metadata.xlsx details the definitions and units for the column headers in the raw data files

README.lichens_as_marine_bioindicators_analysis_procedures.docx details the code and procedures for the data analysis in the paper publication relating to these data (https://doi.org/10.3119/20)