Mesoscale Spatial Patterns of Gulf of Maine Rocky Intertidal Communities

Abstract

Community similarity among macroinvertebrate species assemblages from 12 exposed rocky headlands surveyed in 2004, 2007, and 2012 was examined to resolve mesoscale patterns along an east–west linear distance of 366 km in the coastal Gulf of Maine. The goals were: (1) detect latitudinal patterns of species assemblage similarity and (2) relate species assemblage similarities to environmental factors. Assemblage similarities were correlated with latitude. There was a distinguishable grouping of sampling sites fitting two Gulf regions that separate at mid-coast Maine. This pattern was uniquely intertidal and not shown by subtidal species assemblages. β diversity was high, did not differ between regions, and species turnover accounted for 91% of it. Molluscs and crustaceans, major components of surveyed communities, contributed most of the dissimilarity between regions. Satellite-derived shore and sea surface temperatures explained a significant amount of the variation responsible for producing regional patterns. The regions corresponded with the two principal branches of the Gulf of Maine Coastal Current. These hydrographic features and associated environmental conditions are hypothesized to influence community dynamics and shape the dissimilarity between Gulf regions. The predicted warming of the Gulf of Maine portend change in species turnover from species invasions and range shifts potentially altering rocky intertidal community patterns.

1. Introduction

The delineation of broad scale spatial biodiversity patterns is valuable for detecting, gauging, and predicting the response of communities to environmental change. If variation in community structure produces detectable new patterns, such change can signal modification of community composition with altered or novel species interactions as a consequence [1]. The extent that communities are buffered against change depends on the stability of their populations to recover from perturbations stemming from both environmental and biological factors [2], which in turn will determine the degree of local species extinctions and long-term consequences for community dynamics [3,4,5]. The number and types of potential species interactions in novel communities that emerge can change ecosystem function and linked ecosystem services [6]. Outcomes can have direct economic, demographic, and social consequences for coastal communities especially when commercially valuable species are lost.

The diversity and structure of intertidal communities living at the land–sea interface are shaped by the aggregative effects of local and broad scale marine, terrestrial, and atmospheric processes. What shapes broad scale diversity and complexity of intertidal community structure is intimately tied to coastal circulation and oceanic processes [7,8] and the biogeographic patterns that result are strongly associated with these features [9,10]. The dispersal and delivery of nutrients, food, and propagules are steered by ocean currents, which also set physical limits that intertidal species tolerate. Intertidal communities are also vulnerable to the effects of aerial exposure during periods of low tide. Air temperature, humidity, and precipitation, influenced by broad scale meteorological systems [11], exert selective influences on community composition according to the physiological requirements of species [12]. In addition, many intertidal communities are subject to the near- and far-field effects of freshwater outflow from associated watersheds [13,14,15]. Overall, the shoreline frames a habitat subject to extremes where the effects of climate change from altered terrestrial, atmospheric, and marine environments are concentrated.

Global warming will directly affect thermally sensitive processes, and temperature is a pervasive force on all biological phenomena [16]. Global scale distribution patterns of rocky shore intertidal communities are strongly related to temperature [17]. Latitudinal species distributions are limited by the effects of temperature [18,19], which constrain rates of reproduction and pelagic larval development [20,21]. The thermal challenges encountered by species dispersing into the extremes of their geographic ranges tests their physiological adaptations to the intertidal environment [22,23]. Shifts in species distributional ranges from global warming are limited by the genetic capacity for evolving thermal and phenological adaptations [6,24]. Thus, range shifting may not be continuous progressions in space and time but instead can be punctuated over short time and spatial scales. Local extinctions can happen when intertidal thermal environments surpass the capacity of species to acclimate [25,26] due to limited phenotypic plasticity [12]. Understanding the outcomes of warming is not straightforward since the degree to which changes in temperature effect species interactions is not well understood [25,26]. Species establishment is dependent on a variety of abiotic and biotic factors, such as oceanographic conditions, food limitation, competition, and predation. Given the complexity of abiotic and biotic effects, determining regional patterns in community structure will help to facilitate the prediction of changes from global warming.

The association of coastal circulation and temperature was explored to develop testable predictions concerning the role of these features in structuring intertidal communities within the Gulf of Maine (GoM) (Figure 1). The GoM is described by the expanse of water between Cape Cod, Massachusetts, and southwestern Nova Scotia. Isolated by Browns Bank and Georges Bank from the open Northwest Atlantic, the GoM is a semi-enclosed marginal sea with distinct oceanographic and meteorological features [27,28,29]. The Gulf of Maine Coastal Current (GMCC) is one such feature and is a major influence on the Gulf’s biological productivity [30]. The GMCC receives water from the Scotian Shelf as it flows cyclonically near the 100 m isobath from the Grand Banks to Massachusetts Bay south-west. The GMCC has two principal branches, the Eastern Maine Coastal Current (EMCC), where there is an offshore component, and Western Maine Coastal Current (WMCC). The EMCC extends along the eastern Gulf and flows southwest towards Penobscot Bay and the WMCC originates immediately south of the bay and flows into the southern Gulf [31]. Among other physical characteristics, sea water temperature distinguishes these two currents, the EMCC being colder [32,33]. The general cyclonic flow pattern of the GMCC changes seasonally. During spring, summer, and fall the EMCC flows uninterrupted to mid-coast where upon encountering the Penobscot Bay region a portion flows cyclonically away from the coastline. A portion continues past the mouth of the bay to join the WMCC [31,33]. Flow to the WMCC is regulated by complex hydrographic processes and can range from continuous to complete disruption during different years. Keafer et al. [34] described another GoM hydrographic feature, the low salinity Gulf of Maine Coastal Plume (GOMCP), which lies sandwiched between the coast and GMCC and receives water from major rivers along its southwest flowing course extending from eastern Maine. Winter circulation in the GoM is less well known. However, the cyclonic pattern is less organized and slows down [35,36,37].

The main purpose of this investigation was to examine the similarity of species assemblages among rocky intertidal communities in the GoM to reveal mesoscale spatial patterns and their persistence in time. The strong thermal gradient established by summer coastal circulation was predicted to influence the similarity among species assemblages on northern and southern Gulf shores. This was evaluated using a combination of multivariate and nonparametric approaches to compare patterns in temperatures and community variability across space. On the basis of these analyses, it was determined whether intertidal communities were similar throughout the GoM or if there were regional differences. This led to comparisons with subtidal communities across the same spatial extent. By doing so, differences found between these habitats were used to delineate where community similarities were found.

2. Materials and Methods

2.1. Study Sites and Rocky Intertidal Surveys

The study area spanned approximately 2 degrees of latitude, a distance of 336 km. A total of 12 exposed headlands were surveyed in 2004, 2005, and 2012. Rocky, exposed headlands were selected as study locations to have some degree of habitat similarity for comparisons and give a mesoscale geographic representation of the GoM coast. They were distributed from Sea Point, near the New Hampshire–Maine state border, to West Quoddy Head, near the Maine–New Brunswick, Canada border (Figure 1). In order to keep some degree of congruence among habitats, Central Maine, primarily occupied by Penobscot Bay, was not surveyed because the majority of the exposed rocky shores in the Penobscot Bay estuary are on islands and not the mainland. Estuarine mud covers most of bottom of this island bay complex which receives freshwater from the Penobscot River and its watershed [38]. Intertidal communities were sampled at low tide in summer (June–August) during 2004 with line-transects, and in 2005 and 2012 with walk-about surveys. Time of low water and tidal amplitudes relative to mean lower low water (MLLW) were taken from the WWW Tide and Current Predictor [39]. The presence of macroinvertebrates (≥1 mm) was recorded, identified in the field to the lowest taxon possible, usually species, or if unknown, collected, and identified in the laboratory the same day.

Sampling Methods

Headlands were surveyed in 2004 with line transects extending the full intertidal range from low water (chart datum) occurring at the predicted time to the high water line marked prior to low tide. There were three line transects of equal length for each location and the positions of endpoints recorded with WAAS GPS. Reconnaissance surveys were conducted as part of a pre-selection process for positioning sample transects that best avoided tidepools, large boulders, and upturned bedrock benches. Tide pools were not sampled and when encountered, the meter interval free of standing water closest to the immersed transect sample was selected instead. All macroinvertebrates in every meter interval which contacted a transect line were recorded.

During 2004 surveys, substrate types and dominant algae were assessed in four, non-random 1 m² quadrats positioned along one transect randomly selected from the three line transects. One quadrat was located about 2 m above the lowest exposed point on the shore, another approximately 2 m below high water, and two situated at quarter marks between these stations so that adjacent pairs were equidistant from each other. Within each quadrat, the primary substrate type was identified by visual estimate after dividing each 1 m² quadrat into 0.25 m² subsamples. The substrate class that covered > 50% of the surface was classified as primary. Substrate classes were gravel, cobble, boulder, and rock as defined by Brown [40].

Headlands were surveyed in 2007 and 2012 using walk-about surveys. The area surveyed varied among locations because of differences in slope, shoreline contour, and topography which determined the amount of exposed shore in addition to tidal amplitude. In 2007, each location was sampled over the course of 2 or 3 days during one low tide for 4 h each day. In 2012, each location was sampled in one day during one low tide. The difference in times for completing surveys reflects a time gaining experience with each intertidal site and funding objectives. The procedure for walk about surveys was as follows. Intertidal macroinvertebrates were sampled at randomly selected points with 10 × 10 cm quadrats. These sample points were at the terminus of path segments of random length. Path segments were oriented in randomly chosen compass bearings from a sample point. In this way, each quadrat was an independent, randomly selected sample. Sampling began towards a seaward horizon away from the high tide mark. Upon reaching the water’s edge, the general heading switched to a landward horizon until reaching the high water mark, when the general heading switched back to seaward. Sampling continued until species accumulation curves reached an asymptote. Start and end points of sample paths were recorded with WAAS GPS, landmarks, and photographically. Boundaries to most survey areas were taken from maps accessed from the Critical Areas files in the Maine State Archives Library, otherwise they were defined using ArcMap™.

2.2. Exposure Index

Exposure was estimated for each survey location using an index that combines wind energy and effective fetch [41]. For this index, wind energy (W) depends on the duration and average speed (knots) the wind blows in each compass direction defined by 22.5° sectors. It was calculated using the equation:

$$W=\left(\frac{\mathrm{percentage} \mathrm{of} \mathrm{time} \mathrm{the} \mathrm{wind} \mathrm{blows} \mathrm{in} \mathrm{a} 22.5° \mathrm{compass} \sec \mathrm{tor}}{100}\right)\times {\left(\mathrm{mean} \mathrm{wind} \mathrm{speed}\right)}^2$$

Effective fetch introduces a bathymetric component to the exposure index and is the quotient of actual fetch (F) divided by the sum of the extent (nautical miles) of shallow water <6 m deep joining the shoreline (CS) plus shallow water <6 m deep beyond that margin (DS). Fetch has a 100 NM maximum, i.e., distances greater than 100 NM are recorded as 100 NM. In summary, the exposure index is the sum of wind energy and effective fetch within each 22.5° compass sector of shoreline calculated using the Equation E1 of Thomas [41]:

$${\displaystyle \sum }\log W\times \log \left[1+F/\left(CS+0.1 DS\right)\right]$$

The measurement of each variable was achieved using the following method. Wind roses with 22.5° sectors were generated with WRPLOT View^TM (version 8.0.2) using wind velocity data recorded over a 5-year period between 2004 and 2012 at nearby coastal weather stations. From these, wind duration and mean speed were used to calculate wind energy for each sector and subsequently summed to calculate W. Using Google Earth Pro, wind roses were digitally centered on top of survey locations so that the first compass sector of the rose aligned with true north. After adjusting the transparency of the wind rose, the maximum extents of shallow water (CS and DS) within each sector were measured from NOAA Office of Coast Survey raster navigational charts overlaid on Google Earth imagery.

2.3. Subtidal Species Assemblages

Subtidal species assemblages among GoM benthic communities were explored for patterns of species assemblage similarity to compare with those of intertidal assemblages and the biogeographic analysis by Hale [42]. To carry this out, Environmental Protection Agency National Coastal Condition Assessment (NCCA) [43] data collected from subtidal stations July–September during 2000 to 2004 were selected from the same set of data analyzed by Hale [42] for western Atlantic biogeographic patterns. Station data includes benthic macroinvertebrate species abundance from 0.05 m² grabs, one grab sample per station, with no resampling, and water quality and temperature measurements. Sample mean depth, after removing rivers and ponds, was 19 m (max = 77.9 m, min = 1.1 m, mode = 18 m). Proxy stations were selected for nearness to intertidal survey locations (within 1km) and, when possible, shallow (<10 m) depths. Species presence data were used to characterize subtidal species assemblages. Comparisons with rocky subtidal epifauna were not possible because no data were available for the complete set of intertidal locations and the data that were accessible were collected outside of the time frame of intertidal surveys.

2.4. Coastal Temperature Data Acquisition and Analysis

Coastal land and sea surface temperatures during intertidal surveys were estimated using Copernicus Climate Change Service (C3S) Climate Data Store data. Land temperatures (2 m temperature, i.e., air temperature 2 m above the ground) were downloaded from ERA5—Land monthly data from 1950 to present [44]. The data are monthly averages calculated from average daily temperatures and are gridded with a horizontal resolution of 0.1° × 0.1°. Sea surface temperatures were downloaded from Sea surface temperature daily data from 1981 to present, derived from satellite observations [45]. The chosen Level 4 processing (Version 2.1) yielded temperatures resulting from a combination of measurements made by multiple sensor types (AVHRR, ATSR, SLSTR, and MetOp) and satellites (NOAA, ERS, Envisat, and Sentinel). Gridded SST data have a horizontal resolution of 0.05° × 0.05°. Preliminary examination showed that for all intertidal survey years, the warmest temperatures occurred during July–September and coldest during December–February. Therefore, mean temperatures for the two periods, called summer and winter from hereon, were calculated using contiguous months. For winter, the December of the year preceding January and February was used. For example, 2004 winter temperatures were assembled from 2003 December and 2004 January and February temperatures. Radiometer-based SST temperatures were ground truthed with temperatures measured during surveys and recordings by buoys of the Northeastern Regional Association of Coastal and Ocean Observing System (NERACOOS).

A three-way ANOVA (Sigma Plot 14.5) was used to explore differences in temperatures among survey years, the region where surveys were conducted, and where the temperature was estimated (land versus seas surface) as a factor. Data passed the Shapiro-Wilk and Brown-Forsythe tests for normality and equal variance, respectively. When significance was detected, the Holm-Sidak test was used for multiple comparisons among means to find which were statistically different.

2.5. Statistical Analysis of Species Assemblages

The similarity of species assemblages was compared among survey locations using Plymouth Routines in Multivariate Ecological Research (PRIMER 7) [46,47], PERMANOVA+ [48], and their various subroutines. Only species incidence data were analyzed. Spatial analysis of species distributions within and among locations was not explored. Datasets collected using line transects, walk-about surveys, and NCCA benthic grabs were analyzed separately to accommodate differences in temperatures, surveyed locations, sampling protocols, and year sampled. Species accumulation curves for 2004, 2007, and 2012 surveys showed that all assemblages were adequately sampled with species richness reaching an asymptote (Supplementary Materials, Table S1, Figure S1). Samples were pooled for species presence at each intertidal survey location prior to analysis. Species presence was compiled from species abundance for subtidal NCCA grab samples. Patterns of species assemblage similarity among these four sets of data were investigated using hierarchical cluster analysis, canonical analysis of principle components (CAP), nonparametric multidimensional scaling (nMDS), and tests of mesoscale differences between species assemblages north and south of the mid-coast Penobscot Bay region (ANOSIM). Species similarity within regions and dissimilarity between regions was computed and compared (SIMPER). β diversity and its components was assessed for GoM species assemblages using R. Statistical significance for all tests was defined by p values less than 0.05.

Regional patterns in species assemblage similarity were investigated with hierarchal cluster analysis using the group average as the cluster mode on Bray-Curtis similarity matrices of species presence data. Evidence of statistically distinct clusters was explored with the similarity of profiles test (SIMPROF). An association of species assemblage similarity with latitude was evaluated using the canonical analysis of principle components (CAP). This method was used to visualize the distances between centroids of survey location similarity using latitude as the predictor variable. CAP also assessed the strength of correlation (${\delta }_x^2$) of the constrained ordination of samples with latitude. Patterns in similarity among assemblages were visualized with nonparametric multidimensional scaling (nMDS) on Bray-Curtis similarities. A spatial relationship of assemblage similarity at a coarser scale than latitude was explored by grouping survey locations by region and performing a one-way analysis of similarity (ANOSIM) test, with regions as unordered groups, to evaluate a statistical difference between regions north and south of mid-coast Maine. There were only three replicates, i.e., locations, per group for the 2007 walk-about surveyed headlands, too few to give meaningful ANOSIM results [47]. Instead, a one-way PERMANOVA was used to test for difference with Bray-Curtis similarities and region as a fixed factor in the model, followed by a pair-wise test to resolve statistical differences between north and south regions using Monte Carlo p values (p_(MC)). Unlike ANOSIM, PERMANOVA permutes similarity values rather than ranks, and evaluates the difference between centroids. The problem of a small number of replicates was surmounted by using Monte Carlo p values. Average similarity within and dissimilarity between regions were measured using the similarity percentages routine (SIMPER). This test also named species that contributed most (up to 70%) to the within-region similarity and differences between regions.

β diversity was evaluated for each survey year and for species incidence pooled among years at the Gulf scale using the Sørensen dissimilarity coefficient. Sørensen was chosen since the Bray-Curtis coefficient, used throughout analyses, is identical when calculated on presence/absence data [46]. The contributions of nestedness and turnover to structuring β diversity was assessed by partitioning β diversity into the components of species richness difference and species replacement. These computations were performed using the function beta.multi in the betapart package [49] in R (Version 4.2.0). Next, differences in β diversity between the north and south Gulf regions were explored using the betadisper function in vegan (version 2.6-2) [50]. This analysis used PERMDISP [51] to test if β diversity differed significantly between regions.

2.6. Temperature and Species Assemblage Similarities

The relationships of summer land temperature, summer SST, and exposure with species similarities among surveyed locations were examined. Environmental variables were not strongly collinear (Pearson $\left|r\right|$ ≥ 0.95) and were normalized prior to analyses to place them on a common scale. Summer land and sea surface temperatures and exposure were fitted with Bray–Curtis similarity matrixes using the distance-based linear models (DISTLM) routine in PERMANOVA + [48]. This procedure modelled the relationship between species assemblage similarities using the environmental variables as predictor variables. In general, DISTLM partitions the variation in multivariate data described by a resemblance matrix, and predictor variables are fit individually or sequentially to the model. Thus, the proportion of variation explained by each variable alone (marginal tests) and the proportion explained by each variable when added sequentially to a specified set of variables (conditional tests) are calculated with associated p-values acquired by permutation methods. The conditional sequential tests can disentangle the proportion of variation explained by each variable when added after ones previously fitted to the model. Finally, fitted models were visualized using the distance-based redundancy analysis (dbRDA) routine in PERMANOVA + and the patterns of sample ordination seen on plots examined.

3. Results

3.1. Coastal Temperatures

The thermogeography of the region features temperatures, which vary according to season and latitude. Summer coastal land temperatures were warmer than sea surface temperatures and cooler in the northern Gulf compared to the south (Figure 2). Sea surface temperatures followed this same latitudinal trend. Winter featured coastal land tempera-tures colder than sea surface temperatures (Figure 3). The summer trend with latitude was not present. Year, region, and temperature type (land versus SST) were statistically significant by the three-way ANOVA test, with a significant interaction of region and temperature type (

Table 1. Summary table for three-way ANOVA test of summer land temperatures (2 m air) and SST during intertidal survey years. Region refers to northern versus southern GoM, defined by survey locations. Temperature Type refers to land (2m air) temperature versus sea surface temperature (SST).

Source of Variation	DF		SS	MS	F	p
Year	2		61.256	30.628	40.958	<0.001
Location	1		110.079	110.079	147.207	<0.001
Position	1		187.953	187.953	251.346	<0.001
Year × Region	2		1.41	0.705	0.943	0.398
Year × Temperature Type	2		3.342	1.671	2.234	0.12
Region × Temperature Type	1		20.165	20.165	26.966	<0.001
Year × Region × Temperature Type	2		0.0917	0.0459	0.0613	0.941
Residual	41		30.659	0.748
Total	52		446.112	8.579
Holm–Sidak Pairwise Multiple Comparison Tests
Comparison	Diff of Means		t	p
2012 vs. 2004		2.431	8.888	<0.001
2012 vs. 2007		1.65	5.268	<0.001
2007 vs. 2004		0.78	2.452	0.019
South vs. North		2.995	12.133	<0.001
Land Temperature vs. SST		3.913	15.854	<0.001
Region × Temperature Type Comparison Tests
Comparison	Diff of Means		t	p
Temperature Type within South
Land Temperature vs. SST		2.632	7.294	<0.001
Temperature Type within North
Land Temperature vs. SST		5.195	15.417	<0.001
Region within Land Temperature
South vs. North		1.713	4.88	<0.001
Region within SST
South vs. North		4.277	12.321	<0.001

). Multiple pairwise comparisons found statistical significance: (1) among all survey years, (2) northern and southern regions, and (3) land versus sea surface temperatures. Pair-wise comparisons exploring the interaction between region and temperature type found significance in all combinations of these two factors. In other words, land and sea surface temperatures differed significantly in the north and south, as the north and south regions differed in land temperature and SST.

3.2. Intertidal Description

All of the 12 exposed headlands were bedrock, and many were covered with boulder and cobble in varying degrees (

Table 2. Aspect, exposure, and primary substrate for exposed headland intertidal locations surveyed. Abbreviations for substrate: B, bedrock; Bo, boulder; C, cobble. Symbols show year of survey: *, 2004; ^†, 2007; ^§, 2012.

Survey Site	Location	Aspect	Exposure	Substrate
Sea Point *,†,§	43.09°–70.66°	120° SE	18.48	B/Bo/C
Biddeford Pool §	43.45°–70.33°	116° ESE	75.26	B/Bo
Bailey Island *,†,§	43.72°–70.00°	218° SW	27.98	B
Bald Head §	43.70°–69.85°	245° WSW	62.61	B
Pemaquid Point *,†,§	43.83°–69.51°	190° S	23.69	B
Marshall Point *,§	43.92°–69.26°	122° SE	24.58	B/Bo
Schoodic Point †,§	44.35°–68.08°	230° SW	20.26	B/C
Petit Manan Point *,§	44.40°–67.90°	256° WSW	31.78	B/Bo
Red Head *,†,§	44.45°–67.58°	231° SW	25.18	B
Roque Bluffs *	44.68°–67.15°	198° SSW	24.79	B/Bo
Black Point §	44.68°–67.15°	63° ENE	27.79	B
West Quoddy Head *,†,§	44.81°–66.95°	90° E	32.17	B/Bo

). The exposure index ranged most often from 20–32, although the full range was 75.26–18.48 (median = 26.49). Generally, the dominant macroalgae present were: Fucus vesiculosus, Chondrus crispus, and Mastocarpus stellatus, with Fucus distichus, Saccharina latissima, Laminaria digitata, and Alaria esculenta lowest intertidally. Rocky surfaces were coated in patches with Hildenbrandia sp. and Ceramium sp. mid-intertidally. Cladophora sp. and Corallina officinalis were part of this understory.

3.3. Species Diversity

A pooled total of 117 taxa (Supplementary Materials, Table S2) was dominated by molluscs (29%) and crustaceans (17%). β diversity of species assemblages was moderate among survey years (

Table 3. β diversity and contribution of its components for GoM exposed headland rocky intertidal species assemblages for each survey year and all years pooled.

Year	β diversity	Turnover	Proportion of β	Nestedness	Proportion of β
2004	0.63	0.56	0.89	0.07	0.11
2007	0.49	0.38	0.78	0.11	0.22
2012	0.67	0.59	0.88	0.08	0.09
Pooled	0.86	0.78	0.91	0.08	0.09

). There was no statistical difference in β diversity between south and north regions. Species turnover accounted for 78% to 88% of the β diversity. This indicates that variation among species assemblages results from species replacement along the longitudinal gradient of the GoM shore and not because locations are nested subsets. These trends in β diversity and its components were consistent when species incidence was pooled among years except β diversity was high and species turnover was greater.

3.4. Regional Comparison of Similarity among Exposed Rocky Headland Species Assemblages

Exposed headland species assemblages differed in similarity on a regional scale. For all survey years, assemblages clustered into two statistically distinct groups corresponding to regions north and south of mid-coastal Penobscot Bay. However, there were no significant differences among assemblages within each region (Figure 4). Survey locations based on assemblage similarity clustered by region and latitude in constrained CAP ordinations of species assemblages (Figure 5). The relationship was strong (${\delta }_1^2$, range 0.79–0.94) and canonical correlations were highly significant (

Table 4. Summary of results from the statistical comparison of southern versus northern GoM species assemblage similarities (ANOSIM) and the relationship of species assemblage similarity with latitude (CAP).

Year	ANOSIM	p-Value	$\mathbf{CAP} \left({\mathit{\delta }}_1^2\right)$	p-Value
2004	0.839	0.003	0.93879	0.002
2007	2.12 †	0.03 ††	0.79364	0.1099
2012	0.768	0.002	0.79643	0.002

^† Test statistic, t, for 2007 was calculated from PERMANOVA, not ANOSIM. ^†† Monte Carlo adjusted p.

), except for the 2007 survey due to the low number of locations. Regional dissimilarity was clear for all survey years in unconstrained two-dimensional nMDS ordinations, each with low stress value and resembling CAP ordinations (Figure 6). Southern locations grouped together and separate from northern ones that grouped together. Southern and northern regions differed significantly when assemblage similarities were compared (Table 4). Among all years, the average Bray-Curtis similarity of species present within southern and northern regions ranged among years from 67.21 to 77.3 (

Table 5. Species that contributed most to dissimilarity (up to 70%) between regions determined by SIMPER. Average Bray-Curtis dissimilarity between regions in parentheses. Average Bray-Curtis similarity within regions in braces.

2004		2007		2012
(41.28)		(32.97)		(38.91)
South	North	South	North	South	North
{67.85}	{77.33}	{78.21}	{76.4}	{73.71}	{67.21}
Halichondia panacea	Lacuna vincta	Botryloides violaceus	Scycon ciliarus	Botryloides violaceus	Cucumaria frondosa
Halisarca sp.	Metridium senile	Harmothoe imbricata	Leptasterias littoralis	Botryllus schlosseri	Buccinum undatum
Pectinaria gouldii	Stenosemus albus	Homarus americanus	Aulactinia stella	Didemnum vexillum
Lineus viridis	Urticina felina	Idotea phosphorea		Hemigrapsus sanguineus
	Jaera albifronds	Nereis pelagica		Crepidula fornicata
	Halocynthia pyriformis	Clava multicoirnis		Diplosoma listerium
		Tubularia indivisa		Tubularia indivisa
		Cadlina laevis		Disporella hispida
		Crepidula fornicata
		Hemigrapsus sanguineus
		Onchidoris bilamellata
		Phoxichilidium fematorum

). Average Bray-Curtis dissimilarity between regions ranged from 32.97 to 41.28. Overall, most species which contributed up to 70% of the average dissimilarity were arthropods and molluscs. However, when only the species found exclusively in one or the other region were considered, the dominant taxa changed to a mixed group. In general, more species were found only in the southern region with invasive species appearing in 2007 and 2012. Interestingly, subtidal species assemblages did not cluster according to similarity by region, were not significantly correlated in CAP analysis with latitude (${\delta }_1^2$ = 0.287, p = 0.219), and there was no statistical difference between regional groupings (ANOSIM R = 0.136, p = 0.12). nMDS ordination showed no clear pattern of separation among species assemblages according to where grab samples were taken in respect to regions north and south (Figure 7).

3.5. Spatial Relationship of Communities with Temperature and Exposure

Summer sea surface temperature and land temperature explained a large, statistically significant amount of the variation in intertidal species assemble similarity. When each variable was considered individually in marginal tests, SST explained 37–53% of the variation in species assemblage similarity among all survey years. Land temperatures explained 31–44%. These relationships were statistically significant in all cases except 2007 where only SST was significant (

Table 6. Marginal test results of distance-based linear modelling of species similarity with temperatures and exposure for survey years. Results show the proportion of variability (Prop.) explained for each variable and its level of significance.

Year	Variable	Pseudo-F	p	Prop.
2004	Summer Mean SST	4.873	0.003	0.45
	Summer Mean Land Temp	4.147	0.004	0.41
	Exposure	0.757	0.590	0.11
2007	Summer Mean SST	4.564	0.020	0.53
	Summer Mean Land Temp	3.217	0.062	0.45
	Exposure	0.188	0.922	0.04
2012	Summer Mean SST	5.266	0.001	0.37
	Summer Mean Land Temp	4.076	0.001	0.31
	Exposure	0.783	0.614	0.08

). Exposure did not explain a significant amount of variation in marginal tests (≤11%, p > 0.05). In sequential conditional tests where the ordering of SST and land temperature was switched, the variable explaining the largest proportion of variation was the first one in the sequence evaluated (

Table 7. Sequential conditional test results of distance-based linear modelling of species similarity with temperatures and exposure for survey years. Variables are fit to models as covariables and evaluated as sets. The sequence of SST and land temperature was switched in tests for each year. The proportion of variation (Prop.) of an added variable is its contribution to the amount already explained by the preceding variable (Res. df, residual degrees of freedom).

Year	Variable Test Order	Pseudo-F	p	Prop.	Res. df
2004	+Summer Mean SST	4.873	0.003	0.45	6
	+Summer Mean Land Temp	0.608	0.778	0.06	5
	+Exposure	0.735	0.624	0.08	4
	+Summer Mean Land Temp	4.147	0.001	0.41	6
	+Summer Mean SST	1.010	0.453	0.10	5
	+Exposure	0.735	0.614	0.08	4
2007	+Summer Mean SST	4.564	0.018	0.53	4
	+Summer Mean Land Temp	1.894	0.217	0.18	3
	+Exposure	0.492	0.683	0.06	2
	+Summer Mean Land Temp	3.217	0.070	0.45	4
	+Summer Mean SST	2.808	0.068	0.27	3
	+Exposure	0.492	0.662	0.06	2
2012	+Summer Mean SST	5.266	0.001	0.37	9
	+Summer Mean Land Temp	1.635	0.113	0.11	8
	+Exposure	0.461	0.886	0.03	7
	+Summer Mean Land Temp	4.076	0.001	0.31	9
	+Summer Mean SST	2.512	0.012	0.16	8
	+Exposure	0.461	0.873	0.03	7

). When SST was first, it explained a statistically significant proportion (37–53%) of variation. Land temperatures explained 5–18% more, statistically insignificant amounts. When land temperature was first, it explained more of the variation (31–45%) in species assemblage similarity than SST, statistically significant amounts except for 2007. Sea surface temperature contributed an additional 10–27%, insignificant amounts except for 2012. Placing exposure first in the testing sequence did not change the outcomes for temperatures, and the amount of variation it explained was always the smallest and insignificant. In summary, land temperatures and SST together explained a significant amount of the variability in assemblage similarity but not exposure. The models performed well and captured most of the variation in species assemblage variation as shown by their associated plots produced by distance-based redundancy analysis (Figure 8). Among all years, the first two dbRDA axes explained 90–100% of the fitted variation, which was about 48–78% of the total variation in species assemblage similarities. The separation of species assemblages into northern and southern groups by dbRDA was clear in all plots and consistent among all survey years.

4. Discussion

Exposed rocky headland intertidal species assemblages of northern and southern GoM shores were dissimilar. Penobscot Bay located at mid-coast Maine, the largest estuary in Maine and the second largest on the US east coast [52], marked the division between the two regions. This pattern was persistent among three sets of data, which varied in collection methods and spanned a total of eight years. Variation in species assemblage similarity was correlated with latitude and persistent among surveys separated in time by 8 years. The regional differences were consistent with prior exposed rocky shore studies conducted in similar if not the same locations [53,54]. The dissimilarity was confined to the intertidal since it was not evident among subtidal species assemblages. In general, the northern and southern Gulf regions matched the northern and southern faunal zones proposed by Bousfield and Thomas [55]. Their biogeographic scheme was based on temperature and divides the GoM coastline into three zones: a < 12 °C subarctic zone in northern Maine, a 12–15 °C boreal region in central Maine, and a 15–18 °C cool temperate zone that extends to Massachusetts Bay. Central Maine, primarily occupied by Penobscot Bay, was not surveyed by this study. Adey and Heyek [56] documented differences between northern rocky intertidal communities from those in the southern GoM, work that contrasted the pre-existing idea that the Gulf was part of a single biogeographic unit extending from Cape Cod to Labrador [57]. The pattern of dissimilarity between Gulf regions is not limited to exposed rocky shores. It was found among macroinvertebrate species assemblages of sand beaches [58] and low energy intertidal areas dominated by Ascophyllum nodosum [59]. The pattern of dissimilarity was a key feature of GoM intertidal communities with a signal strong enough to be detected by species presence data.

Species were not evenly distributed among all surveyed locations, which lead to the dissimilarity between north and south regions. Species turnover contributed most to this pattern. Analysis of these dissimilarities with SIMPER showed trends that coincided with GoM invasive species histories, each with sea water temperature identified as a common factor underlying population dynamics. Some species were found exclusively in one region, a situation that occurred more often in the southern Gulf, particularly in 2007 and 2012 when invasive species accounted for some of them. These included the colonial tunicates Botryloides violaceus, Didemnum vexillum, and Diplosoma listereum, and Botryllus schlosseri currently understood to be cryptogenic [60]. Around the 2007, the southern region was subject to invasion by those species and their community dynamics, competition for space in particular, were shown to be correlated with seasonal changes in sea water temperatures [61]. Similarly identified among the southern Gulf exclusives was the invasive Asian shore crab Hemigrapsus sanguineus that appeared in the region around 2001 [62]. Based on 2002 to 2005 coastal Maine surveys that over-lapped the range in latitude of the current study, Stephenson et al. [62] suggested that sea water temperatures might limit range expansion to the warmer coasts south of Penobscot Bay. Since then, H. sanguineus has increased in density in the southern region where established populations remain confined [63].

Shore and sea surface temperatures explained a large proportion of the variation of similarity among species assemblages. When evaluated individually or as covariates, the degree of variation explained was statistically significant except when the number of surveyed locations was small, i.e., the 2007 surveys. The influence of shore temperatures reinforced the idea that regional community dissimilarity is a feature of the shore and not the subtidal. Their contribution was significant despite the lower resolution of the gridded data compared to SST. The different thermogeography of northern and southern Gulf regions is largely influenced by ocean circulation, and the flow of major coastal currents match the pattern of dissimilarity among species assemblages of the two regions. Coastal water temperatures do not appear to drive land temperatures [64], but they do act as a buffer, cooling coastal land temperatures in summer and warming them in winter [26,27]. Across-shore thermal gradients that are discontinuous in summer could sort community composition according to species thermal tolerance to produce the patterns in species similarity. Elsewhere, sea surface temperatures were found to play key roles in the regional distribution patterns of species assemblages [9,10,65,66,67]. Likewise, shore temperatures influenced species intertidal distributions [1,12,23,68]. Temperature can affect major ecological patterns and community assembly by driving metabolism, resulting in modifications of longevity, population growth, and consequently competition [16] among other species interactions [25,69]. The effects of temperature are pervasive throughout biological processes, and temperature zonation and biogeographic regions are some macroscale manifestations [18,70].

The lack of a significant effect from wave exposure was unexpected since it is known to influence species richness and diversity among intertidal communities in the same biogeographic region [57,61,62,63,64,65,66,67,68,69,70,71,72,73] and elsewhere [74]. The method for calculating exposure was not an issue since it was supported by survey data [41,75]. This suggests that the effects of exposure are best explored using abundance data and vertically stratified sampling methods to include shore height. In addition, the range in exposure indices was not very large especially when the three extreme measures were ignored (range: 20.26–32.17; SD = 3.86). Therefore, the variation in exposure among survey locations may be too slight to evaluate a relationship with assemblage similarity. Future studies examining an effect of wave exposure on intertidal communities might include sheltered and partially exposed shores [74] to increase the range of variation of exposure among survey locations.

Rocky intertidal communities on northern and southern shores are dissimilar and the mid-coast Penobscot Bay region marks the area where the two regions diverge. However, other research indicates that this feature is not expected. Within the GoM, population genetics of some of the same species found in surveys show no discontinuities [76,77]. Instead, a review of population genetics data [77] showed a discontinuity displaced to the south of Cape Cod, not within the GoM. However, since that review, two species of rocky intertidal gastropods with non-planktonic development, Nucella lapillus [78] and Littorina saxatilis [79], were shown to split into northern and southern clades within the same two regions defined by the current study. Models integrating ocean currents and species with high-dispersal larvae did not predict a peak of range boundaries within the GoM [80]. That said, oceanographic features of the northern Gulf, the EMCC in particular, appear to set the southern range boundary for the mussel Mytilus trossulus [81]. Large-scale analyses of West Atlantic species distribution patterns support conclusions from population genetics and modelling [82]. Furthermore, Hale [43] did not find a transition area for subtidal benthic invertebrates at mid-coast Maine, a finding corroborated by the present study, but instead found one to the south of Casco Bay. That embayment is located approximately 90 km south of the mid-coast. In summary, there is much evidence to the contrary of a mid-coast Maine rocky intertidal discontinuity.

How might the differences between predicted boundaries and the results of the present study be reconciled? Firstly, the transition area for benthic species assemblages was identified for species that are subtidal soft bottom inhabitants, not rocky intertidal ones. Additionally, the species assemblages were different. Subtidal assemblages were dominated by polychaetes, with crustaceans and molluscs the next most abundant taxa [43]. Rocky intertidal assemblages were dominated by molluscs and crustaceans. Next, the incongruity stemming from population genetics has value in the sense that these studies rule out the possibility that hydrodynamic barriers to gene flow via larval dispersal shape regional divergence. Genetic differences among populations are not a prerequisite for community dissimilarity. Species interactions can influence community assembly [83]. Likewise, the meta-analysis of species distributions [82] did not consider community dynamics and its consequences for predicting discontinuities. Finally, high-dispersing larvae are not characteristic of all rocky intertidal species so the hydrodynamic modelling of discontinuities [80] is limited. These models also did not consider how community dynamics might influence their predictions. While species interactions were not specifically examined by the present study, they are implicit in structuring the similarities among the surveyed species assemblages. The dissimilarities between northern and southern shores are likely a signal of the interactions particular to the sets of species present in those regions, an idea supported by their differences in predation [84], recovery from disturbance [54], and possibly recruitment [81,85].

The GoM is currently warming faster than most other bodies of waters globally [86,87]. The survey conducted in 2012 occurred during an ocean heat wave when warming was especially pronounced in the GoM [86]. In the decade since, there have been profound consequences for fisheries [86,88,89], kelp forests [90], and phenologies [91,92,93]. Given the predicted conditions in the GoM [94], change in species turnover from species invasions [95] and range shifts [96] portend novel species interactions with consequences for rocky intertidal community patterns.

5. Conclusions

GoM rocky intertidal communities were similar within regions, but the regions were distinct and located south and north of the Penobscot Bay estuary. This discontinuity did not extend into the subtidal; it was a uniquely intertidal feature. Both land and SST explained a significant amount of the variation which gave rise to regional dissimilarity, but they did not explain all of it. Species interactions and community dynamics are predicted to play important roles.