Biogeographic patterns of community diversity associated with an introduced alga

Non-native foundation species may alter physical environments and provide habitat, thereby impacting recipient communities. Along the US east coast, we assessed biogeographic patterns of free-living and parasitic community diversity associated with the non-native red alga Gracilaria vermiculophylla, which is characterized by fixed (with holdfast) or free-floating thalli depending on the availability of hard substratum. In summer 2019, we surveyed 17 sites across 3 biogeographic regions. We used a random quadrat design to collect G. vermiculophylla and associated mobile macroinvertebrates per site, and we took abiotic measurements. We also haphazardly collected 100 Ilyanassa obsoleta snails per site to assess trematode diversity. In the lab, macroinvertebrates were removed from thalli and identified to lowest taxonomic level, and snails were dissected to determine trematode prevalence and diversity. Biotic and abiotic variables were analyzed for the best sets of predictors for species richness, abundance, and diversity of macroinvertebrates and trematodes across bioregions. Gracilaria vermiculophylla biomass was used as an offset in free-living analyses. Across all our US east coast sites, we detected 10,113 free-living (mobile) macroinvertebrates across 39 taxa. Three Gammaridean amphipods (Gammarus mucronatus, Ampithoe longimana, and Gammarus lawrencianus) comprised >50% of all detected organisms. We found biogeographic region to be a key predictor of macroinvertebrate abundance and richness. Trematode prevalence and richness were best explained by G. vermiculophylla biomass, while biogeographic region best explained diversity. As a widespread invader, our study provides evidence for associations that have formed as this foundation species has become established outside its native range. Over time, the presence and spread of G. vermiculophylla could continue to impact macroinvertebrate structure and diversity, and future work should directly compare macroinvertebrate communities with G. vermiculophylla to other foundation species along coastlines it is now common. 

Methods
Study Sites: We identified sample sites with verified G. vermiculophylla presence from previous studies (Nettleton et al. 2013; Krueger-Hadfield et al. 2017). In summer 2019, we sampled 17 U.S. East Coast sites, capturing much of the species’ introduced range and encompassing two major geographic barriers at Cape Hatteras and Cape Cod (Engle & Summers 1999; Spalding et al. 2007; Hale 2010). Since summer temperatures are lagged in northern versus southern latitudes, southern sites were sampled earlier than northern sites (Table S1): South of Cape Hatteras (henceforth, SCH) was sampled May 21 – June 22 (n=6, avg=29.9oC); Virginian Province (between Cape Hatteras and Cape Cod: henceforth, VP) was sampled June 25 – Aug 3 (n=8, avg=26.7oC); and North of Cape Cod (henceforth NCC) was sampled Aug 4 – 9 (n = 3, avg=23.5oC). Water temperatures were taken in the shallow intertidal zone just before low tide (see below). For this study, since we were specifically interested in conducting a biogeographic study of macroinvertebrates along the temperate coastline of the U.S. east coast, and given the large swath of sites within VP vs. other biogeographic regions as well as the availability of sites with verified G. vermiculophylla presence/accessibility, the number of study sites per biogeographic region was uneven.

Sampling of Associated Free-Living Macroinvertebrates: We sampled each site for G. vermiculophylla in the shallow intertidal zone while thalli were still submerged before low tide. At each site, we established a 30-meter transect tape along the water-land interface, selected five random numbers (1-30) using a random number generator (each number representing a meter marker on the 30-meter transect), and collected all G. vermiculophylla clumps from those five randomly selected 0.25 m2 quadrats along the transect. We sampled environmental parameters (water temperature, salinity) using a handheld YSI Pro-1030 (Yellow Springs, OH).

We placed sealed bags of algae and water immediately into coolers and then transported them to the lab for processing. In the lab, we soaked the G. vermiculophylla from each replicate in a large bin filled with fresh tap water to induce osmotic shock in the associated mobile macroinvertebrates (e.g., Blakeslee et al. 2016; Fowler et al. 2016). We then used a Fisher Scientific™ 250 micron sieve to separate macroinvertebrates from macroalgae; upon separation, we preserved macroinvertebrates in Pharmaco™ 200 proof Ethyl Alcohol. After shaking off excess water, we weighed the thalli to obtain wet weights (g).

Following field surveys at all sites, macroinvertebrates were dyed with Rose Bengal (Gbogbo et al. 2020) and identified to the lowest possible taxonomic level using guidebooks and keys (Bousfield 1973; Johnson & Allen. 2012). Organisms were observed using a Zeiss MS Series Fixed Magnification Stereo Microscope (6x) and/or a Neatfi Elite XL HD Magnifying Lamp (5x). Gammaridean amphipods, which comprised up to 75% of the total macroinvertebrates at sites (see Results), can be difficult to identify to species level using morphology alone. We therefore classified amphipods into morphotypes and then later barcoded those morphotypes using standard DNA protocols (e.g., Blakeslee et al. 2020a). This allowed us to identify amphipods to species level by BLASTing our resultant sequence data for each morphotype using the NCBI database (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

Sampling of Trematode Parasites: We collected all I. obsoleta at the same sites as described above, except for Provincetown, MA, where I. obsoleta was not found (parasite data = 16 sites). We used the same 30-meter transect tape and 0.25 m2 quadrats as above to collect snails; however, G. vermiculophylla and I. obsoleta were placed into separate bags. We counted all I. obsoleta per quadrat, and then randomly selected 100 snails across the five quadrats to dissect. Because birds are common final hosts for trematode parasites, we also counted the total number of birds by species (waders, seabirds, and dabblers) at each site using a point-count method, while standing stationary for 10 minutes (Byers et al. 2008). In the lab, we measured each live gastropod using digital calipers and then dissected gonad tissues under a Zeiss™ MS Series Fixed Magnification Stereo Microscope at 6x magnification. If infected, we identified the digenean trematode to species level based on its rediae/sporocyst and cercarial morphology using published images and keys and prior knowledge within the lab (Curtis & Hurd 1983; Curtis 1985; Esch et al. 2001; Blakeslee et al. 2012). At the snail stage, trematodes asexually produce hundreds to thousands of clones on a continual basis (i.e., the “reproductive firepower” of digenean trematodes; Rohde 2005); thus we did not count cercariae, rediae, or sporocysts within infected snails, as these do not represent genetically-distinct individuals (e.g., Blakeslee & Byers 2008).

Statistical Analyses: To explore which factors best explained the patterns in the communities we observed, we used Generalized Linear Mixed Models (GLMM) in R 4.2.2 (package glmmTMB) for free-living macroinvertebrates (using site as a random effect, families = Gaussian for all response variables) and Generalized Linear Model (GLM) for parasites (families: prevalence = Binomial, richness = Poisson, diversity = Gaussian). For parasite analyses, we used GLM models that included fixed effects only, because at each site, there were no replicates of response variables, since n=100 snails were selected to be dissected randomly across all replicates. Due to the unevenness in detecting fixed versus free-floating G. vermiculophylla across sites, we did not have the number of replicates to analyze algal type (fixed or free-floating) as a predictor in our statistical models; as a result, we separately analyzed the abundance, richness, and diversity of macroinvertebrates associated with fixed and free-floating G. vermiculophylla thalli using two-tailed t-tests across all sites.

For free-living macroinvertebrates, we identified a strong positive and significant relationship between invertebrate raw counts, richness (number of species), and diversity index (Shannon-Wiener Diversity Index) with G. vermiculophylla biomass (Figure 3). To prevent overparameterization and to standardize the response variables, we adjusted our response variables with G. vermiculophylla biomass as an offset [abundance: square root(raw count/G. vermiculophylla biomass); richness: log(number of species/G. vermiculophylla biomass + 1); diversity: square root(Shannon-Wiener Diversity Index)]. We applied these transformations for these response variables to avoid zero variance problems in our models. For parasites, since we dissected an equal number of snails per site (n = 100) and all the three response variables were obtained from those individuals, we did not need to standardize by G. vermiculophylla biomass. For parasites, prevalence was the proportion of infected I. obsoleta out of 100 randomly dissected snails per site; richness was the number of digenean species; and diversity was the Shannon-Weiner Diversity Index.

We selected biologically-relevant predictors for our models after testing for autocorrelations. For free-living organisms, these predictors were water temperature (°C), salinity (PPT), and biogeographic region, with site as a random effect. Since our response variables were standardized by biomass of G. vermiculophylla, the biomass of the seaweed was not included as a predictor in these models to reduce overparameterization. We constructed rarefaction and extrapolation curves to determine the expected number of species per biogeographic region across accumulated individuals using EstimateS (v 9.1.0). For parasites, the predictors were G. vermiculophylla biomass (g), water temperature (°C), salinity (PPT), average snail count, seabird and wading bird count, and biogeographic region.

For both free-living and parasitic communities, we used the corrected Akaike’s Information Criterion (AICc) to determine which model, or sets of environmental variables, best explained the dependent variables of free-living macroinvertebrates and parasites (package AICcmodavg). AICc compares multiple models with different combinations of independent variables (Anderson & Burnham 2002). We used DAICc of ≤2.0 as a cutoff value to determine the top models. Based on our AICc results and their selected predictors in top performing models, we conducted a series of univariate analyses to observe how response variables change for both free-living macroinvertebrates and parasites with key predictors. 

For free-living macroinvertebrates, we used a Nonmetric Multidimensional Scaling (nMDS) to create a two-dimensional ordination plane to visually evaluate community composition among sites (Clarke & Warwick 2001). Per recommendations by Cao et al. (2001), we removed species that occurred <5% in nMDS analyses, and we used square-root transformation and Bray-Curtis Similarity (Clarke & Warwick 2001). For free-living and parasite organisms, we also conducted Similarity of Percentage (SIMPER) analyses to determine the percent each species contributed to the differences observed between biogeographic regions (Clarke 1993; Clarke & Warwick 2001). For free-living macroinvertebrates, we used abundance standardized by G. vermiculophylla biomass. These latter analyses and figures were created using PRIMER-e (v.7).

非本土建群种可能改变物理环境并提供栖息地，进而影响受纳群落。沿美国东海岸，我们评估了与非本土红藻细基江蓠（Gracilaria vermiculophylla）相关的自由生活和寄生群落多样性的生物地理格局——该物种依据硬质底质的可获得性，呈现固着（具固着器）或自由漂浮的藻体。2019年夏季，我们在3个生物地理区域内的17个采样点开展了调查。我们采用随机样方设计，在每个采样点采集细基江蓠及其附着的移动型大型底栖无脊椎动物，并开展非生物环境测量。同时，我们在每个采样点随机采集100枚光滑狭口螺（Ilyanassa obsoleta），以评估吸虫多样性。实验室中，我们从藻体上分离出大型无脊椎动物并鉴定至最低分类阶元，同时解剖螺类以确定吸虫的感染率和多样性。针对生物和非生物变量，我们分析了能够最佳预测跨生物地理区域的大型无脊椎动物和吸虫的物种丰富度、丰度及多样性的变量集。在自由生活类群的分析中，我们以细基江蓠生物量作为偏移量。

在所有美国东海岸采样点中，我们共记录到39个类群的10113只自由生活（移动型）大型无脊椎动物。其中3种端足类钩虾（Gammarus mucronatus、Ampithoe longimana和Gammarus lawrencianus）占所有记录个体的50%以上。我们发现生物地理区域是影响大型无脊椎动物丰度和物种丰富度的关键预测因子。吸虫的感染率和物种丰富度最佳由细基江蓠生物量解释，而生物地理区域则最佳解释了其多样性。作为一种广泛分布的入侵物种，本研究为该建群种在原生分布区外定殖后形成的关联提供了证据。随着时间推移，细基江蓠的存在和扩散可能持续影响大型无脊椎动物的群落结构和多样性，未来研究应直接对比细基江蓠附着的大型无脊椎动物群落与该海岸带其他建群种的相关群落。

## 方法

### 研究样地
我们从既往研究中筛选出经证实存在细基江蓠的采样点（Nettleton et al. 2013; Krueger-Hadfield et al. 2017）。2019年夏季，我们在美国东海岸共采样17个点位，覆盖了该物种的大部分引入分布范围，并囊括了哈特拉斯角和科德角这两个主要地理屏障（Engle & Summers 1999; Spalding et al. 2007; Hale 2010）。由于高纬度与低纬度地区的夏季温度存在滞后性，南部点位的采样早于北部点位（补充表S1）：哈特拉斯角以南（下称SCH）的采样时间为5月21日至6月22日（n=6，平均水温29.9℃）；弗吉尼亚省（哈特拉斯角与科德角之间，下称VP）的采样时间为6月25日至8月3日（n=8，平均水温26.7℃）；科德角以北（下称NCC）的采样时间为8月4日至9日（n=3，平均水温23.5℃）。我们在低潮前于浅潮间带测量了水温。鉴于本研究聚焦美国东海岸温带海岸带大型无脊椎动物的生物地理研究，且弗吉尼亚省点位数量远多于其他生物地理区域，同时结合经证实存在细基江蓠的点位可及性，各生物地理区域的采样点数量并不均等。

### 附着自由生活大型无脊椎动物的采样
我们在低潮前藻体仍处于浸没状态时，于浅潮间带采集细基江蓠。在每个采样点，我们沿水陆交界设置一条30米的样带，使用随机数生成器选取5个1-30之间的随机数（对应样带上的米标位置），并从沿样带的5个随机选取的0.25㎡样方中采集所有细基江蓠藻团。我们使用手持式YSI Pro-1030水质分析仪（Yellow Springs, OH）测量环境参数（水温、盐度）。

将装有藻体和海水的密封袋立即放入冷藏箱，转运至实验室处理。实验室中，我们将每个重复的细基江蓠藻体浸泡在盛有新鲜自来水的大容器中，通过渗透胁迫使附着的移动型大型无脊椎动物脱离藻体（Blakeslee et al. 2016; Fowler et al. 2016）。随后使用Fisher Scientific™ 250μm筛网分离大型无脊椎动物与大型藻类，分离后将无脊椎动物保存在Pharmaco™ 200度无水乙醇中。沥干多余水分后，称量藻体的湿重（单位：g）。

所有点位的野外采样完成后，我们用玫瑰红染料对大型无脊椎动物进行染色（Gbogbo et al. 2020），并通过分类图鉴和检索表鉴定至最低分类阶元（Bousfield 1973; Johnson & Allen 2012）。观察使用Zeiss MS系列固定放大率体视显微镜（6倍）和/或Neatfi Elite XL HD放大灯（5倍）。占采样点总大型无脊椎动物多达75%的端足类钩虾（见结果部分）仅通过形态学难以鉴定至物种水平，因此我们先将钩虾划分为形态型，随后通过标准DNA测序方案对这些形态型进行条形码鉴定（Blakeslee et al. 2020a）。通过将得到的序列数据在NCBI数据库（https://blast.ncbi.nlm.nih.gov/Blast.cgi）中进行BLAST比对，我们可将钩虾鉴定至物种水平。

### 吸虫寄生虫的采样
我们在上述相同点位采集光滑狭口螺，但马萨诸塞州普罗文斯敦点位未发现该螺类，因此寄生虫数据仅来自16个点位。我们使用与前述相同的30米样带和0.25㎡样方采集螺类，但将细基江蓠和光滑狭口螺分别装入不同袋子。统计每个样方内的光滑狭口螺数量，随后从5个样方中随机选取100枚螺类进行解剖。由于鸟类是吸虫的常见终末宿主，我们采用点计数法，在每个点位静止站立10分钟，统计所有鸟类的种类和数量（涉禽、海鸟和钻水鸭）（Byers et al. 2008）。实验室中，我们使用数显游标卡尺测量每只活腹足类的壳长，随后在6倍放大率的Zeiss™ MS系列固定放大率体视显微镜下解剖其性腺组织。若发现感染，我们根据其雷蚴/胞蚴和尾蚴的形态，结合已发表的图谱、检索表以及实验室内部的前期经验，将复殖目吸虫鉴定至物种水平（Curtis & Hurd 1983; Curtis 1985; Esch et al. 2001; Blakeslee et al. 2012）。在螺类宿主体内，吸虫通过无性繁殖持续产生数百至数千个克隆个体（即复殖目吸虫的“繁殖火力”；Rohde 2005），因此我们不对感染螺体内的尾蚴、雷蚴或胞蚴进行计数，因为这些个体并非遗传上独立的个体（Blakeslee & Byers 2008）。

### 统计分析
为探究何种因素能够最佳解释我们观测到的群落格局，我们使用R 4.2.2软件中的广义线性混合模型（Generalized Linear Mixed Models, GLMM）（glmmTMB包）分析自由生活型大型无脊椎动物（以采样点作为随机效应，所有响应变量均采用高斯分布族），并使用广义线性模型（Generalized Linear Model, GLM）分析寄生虫数据（分布族：感染率为二项分布，物种丰富度为泊松分布，多样性为高斯分布）。寄生虫分析仅使用固定效应模型，因为每个采样点仅从所有重复样方中随机选取100枚螺类进行解剖，无重复响应变量。由于各采样点中固着型和自由漂浮型细基江蓠的检出数量不均，我们没有足够的重复样本将藻体类型（固着或自由漂浮）作为预测变量纳入统计模型。因此，我们针对所有采样点，分别对与固着型和自由漂浮型细基江蓠藻体相关的大型无脊椎动物的丰度、物种丰富度和多样性进行双样本t检验。

对于自由生活型大型无脊椎动物，我们发现无脊椎动物的原始计数、物种丰富度（物种数）和香农-威纳多样性指数（Shannon-Wiener Diversity Index）与细基江蓠生物量之间存在显著的正相关关系（图3）。为避免模型过参数化并标准化响应变量，我们以细基江蓠生物量作为偏移量对响应变量进行调整：丰度：√(原始计数/细基江蓠生物量)；物种丰富度：ln(物种数/细基江蓠生物量 + 1)；多样性：√(香农-威纳多样性指数)。我们对这些响应变量进行上述转换以避免模型中出现零方差问题。对于寄生虫数据，由于每个采样点均解剖了等量的螺类（n=100），且三个响应变量均来自这些个体，因此无需通过细基江蓠生物量进行标准化。寄生虫的感染率为每个采样点中100枚随机解剖螺类的感染比例；物种丰富度为复殖目吸虫的物种数；多样性为香农-威纳多样性指数。

我们在检验自相关性后，选择具有生物学意义的预测变量。对于自由生活类群，预测变量包括水温（℃）、盐度（PPT）和生物地理区域，以采样点作为随机效应。由于我们的响应变量已通过细基江蓠生物量进行标准化，因此未将藻体生物量作为预测变量纳入模型，以避免过参数化。我们使用EstimateS（v 9.1.0）构建稀疏和外推曲线，以确定累积个体数下每个生物地理区域的预期物种数。对于寄生虫数据，预测变量包括细基江蓠生物量（g）、水温（℃）、盐度（PPT）、平均螺类数量、海鸟和涉禽数量以及生物地理区域。

对于自由生活和寄生群落，我们均使用校正后的赤池信息准则（Akaike’s Information Criterion, AICc）（AICcmodavg包）来确定能够最佳解释自由生活型大型无脊椎动物和寄生虫响应变量的模型或环境变量组合（Anderson & Burnham 2002）。AICc可对比不同自变量组合的多个模型。我们以ΔAICc ≤2.0作为筛选最优模型的阈值。基于AICc结果和最优模型中的预测变量，我们开展了一系列单变量分析，以观察自由生活型大型无脊椎动物和寄生虫的响应变量如何随关键预测变量变化。

对于自由生活型大型无脊椎动物，我们使用非度量多维标度（Nonmetric Multidimensional Scaling, nMDS）构建二维排序平面，以可视化评估各采样点间的群落组成差异（Clarke & Warwick 2001）。根据Cao等人（2001）的建议，我们移除了在nMDS分析中出现频率低于5%的物种，并使用平方根转换和Bray-Curtis相似性指数（Clarke & Warwick 2001）。对于自由生活类群和寄生虫，我们还开展了相似性百分比（Similarity of Percentage, SIMPER）分析，以确定每个物种对生物地理区域间差异的贡献百分比（Clarke 1993; Clarke & Warwick 2001）。对于自由生活型大型无脊椎动物，我们使用经细基江蓠生物量标准化后的丰度数据。上述分析及图表均通过PRIMER-e（v.7）软件完成。