福寿螺(Pomacea canaliculata)属软体动物门(Mollusca)腹足纲(Gastropoda)中腹足目(Mesogastropoda)瓶螺科(Aillpullaridae), 2000年被世界自然保护联盟(IUCN)列为全球100种最具破坏性的入侵物种之一^[1]。福寿螺原产于南美洲亚马逊流域, 自20世纪80年代被引入我国, 之后广泛扩散, 现最北已扩散至辽宁省朝阳市、最西至甘肃省兰州市^[2—4]。福寿螺具有两栖性、食性广、适应性和繁殖力强等特点, 可迅速建立庞大种群, 对入侵地的农业生产和生态系统完整性等造成严重破坏^[5]。2000年, 广西1万hm²农田受害, 2006年受害面积更达到430万hm², 导致水稻减产10.3万吨, 经济损失1764万元^[6]。此外, 福寿螺还是广州管圆线虫的中间宿主, 一只福寿螺体内可检测出6000条广州管圆线虫, 生食或食用未煮熟的螺肉, 极易引起嗜酸性粒细胞脑膜炎^[7]。长期以来, 相关部门对福寿螺进行大规模防控, 虽在某种程度上抑制福寿螺增长进程, 但效果不显著。全球气候变暖等因素, 导致福寿螺正在以每年8—10 km的速度向北扩散^[8]。

有关福寿螺的分布现状^[9]、繁殖力大小^[10]、环境与福寿螺之间的关系^[11]等方面已开展了较为系统的研究, 这些研究表明福寿螺在人为干扰较大的地区分布较广, 其密度和繁殖力等受多种环境因子的共同影响^{[9, 12]}。然而, 有关福寿螺的种群结构分布特征及其驱动机制的认知明显不足, 主要体现在: (1)研究的空间尺度较小。关于福寿螺的调查大多仅聚焦在单次采样和中小型水体, 鲜有涉及流域尺度上的调查。(2)福寿螺的基本生物学和种群生态学知识的匮乏。(3)未考虑福寿螺不同生活史时期对环境和空间变化的敏感性。

本研究以亚热带桂江流域为例, 基于枯平丰三个时期调查, 从福寿螺种群结构和生活史方面解析其生物学特性, 量化福寿螺种群结构时空分布特征, 并探索福寿螺种群结构与环境、空间因子之间的关系及其驱动机制, 以期为了解桂江流域福寿螺种群现状和预测福寿螺种群未来动态提供研究基础。

1 研究区域与研究方法 1.1 桂江流域概况和样点设置

桂江流域(23°28′—25°55′N, 110°23′—111°34′E)地处亚热带季风气候区, 主源为六洞河, 发源于广西省兴安县境内海拔2000 m的猫儿山。桂江干流全长426 km, 流域面积19288 km²。河流自西北向东南方向流动, 途经兴安县、桂林市、阳朔县、平乐县和昭平县, 于梧州市注入西江干流。上游为源头至桂林市水文站, 地质地貌主要为花岗岩、碎屑岩地层; 中游为桂林市至平乐县, 主要为石灰岩地层; 下游为平乐县至河口, 以砂岩地层为主^[13]。

桂江流域多年平均气温18—21℃, 年平均降雨量1500—2500 mm, 主要集中于3—8月, 占全年降雨量的78%;4—6月降雨量占全年的50%以上^[14]。桂江流域气候温和、雨量充沛、光照充足、水系发达、地貌复杂, 存在较大的环境梯度和空间异质性, 是研究福寿螺的天然实验室^[6]。本研究在桂江流域共设置51个样点(上游18个、中游19个、下游14个)(图 1), 分别于2021年12月、2022年4月和8月进行3次采样。

1.2 样品采集和处理

福寿螺定量采集使用3 m × 3 m的样框^[15], 每个样点采集3次重复样, 将样框中所有福寿螺个体置于封口袋中带回实验室。测量的可量性状包括螺重、壳高、壳宽、壳口宽和壳口长。采用电子秤测量福寿螺体重(精确到0.1 g), 称量前先用滤纸吸干螺体表的水分; 采用量程为20 cm的数显游标卡尺测量壳高、壳宽、壳口宽和壳口长(精确到0.01 mm)(图 2)。性别比例是指种群中雄性与雌性个体数的比例。雌雄可通过观察外部形态(雌螺的厣平滑略凹陷, 雄螺的厣外凸且边缘向内微凹)和内部结构(雄螺存在明显的生殖器)进行鉴定。

1.3 环境因子的测定

现场用YSI便携式水质分析仪(YSI 6680)测定水温、溶氧、pH等环境因子, 用手持式GPS仪(Garmin GPS-76)记录各样点的地理坐标, 用便携式测深仪(Speedtech SM-5A)测量水深, 用流速仪(So LS1206B)测量流速, 用测距仪(Leica YS-TB-CJYONE)测定河宽。水样用1 L采水器采取表层水, 当天带回实验室进行化学参数测量。根据《水和废水监测分析方法(第四版)》^[16], 测定了6个水化学参数, 包括总氮、总磷、氨氮、硝态氮、亚硝态氮、高锰酸盐指数(CODMn)。

1.4 数据分析

用双因素方差分析(two-way ANOVA)检测时间(枯平丰期)和空间(上中下游)对福寿螺种群参数的影响。若组间差异显著, 则进行多重比较, 如果方差同质, 则进行LSD检验, 如果方差不同质, 则选用Games-Howell检验^[17]。为提高数据的正态性和方差齐性, 所有生物数据和环境因子(pH除外)进行lg(x+1)转换。将福寿螺个体按照壳高划分为10组(A—J组): A=0—20.00 mm、B=20.00—24.99 mm、C=25.00—29.99 mm、D=30.00—34.99 mm、E=35.00—39.99 mm、F=40.00—44.99 mm、G=45.00—49.99 mm、H=50.00—54.99 mm、I=55.00—59.99 mm、J>60.00 mm, 壳高小于30 mm的福寿螺个体为幼螺(A—C), 30—50 mm为中螺(D—G), 大于50 mm为成螺(H—J)^[18]。用双因素相似性分析(two-way ANOSIM)检验福寿螺种群结构的时空分布差异。用克里格插值法(Kriging)描述福寿螺时空分布。双因素方差分析、双因素相似性分析和克里格插值法分别在软件SPSS 26.0、PRIMER 6.0和ArcGIS 10.8中完成。

采用主轴邻距法(PCNM)模拟样点间的空间关系, 作为扩散过程的有效替代。PCNM是基于主轴分析的样本间空间距离的结尾矩阵实现的, 选取具正特征值的空间变量, 特征值越小代表的空间尺度越大^[19]。最终, 枯平丰期分别得到了16、18、19个具有正特征值的空间变量。空间因子采用R软件(4.2.1)vegan包中的pcnm函数获得。

采用约束排序分析方法解析桂江流域福寿螺种群结构与驱动因子间的关系。首先, 对福寿螺密度数据进行Hellinger转化, 对环境数据(pH除外)进行lg(x+1)转换, 以保证数据的正态性和方差齐性, 删除相关性较高(|r|＞0.75)和变异波动系数较大(VIF>20)的环境变量。然后, 对福寿螺种群进行去趋势对应分析(Detrended correspondence analysis, DCA), 根据最大轴长选择合适的模型(线性模型(Redundancy analysis, RDA)或单峰模型(Canonical correspondence analysis, CCA)), 利用前向选择(Forward selection)筛选出关键环境因子和空间因子^[20]。最后, 利用变差分解(Variation partitioning)计算出关键环境和空间因子对福寿螺种群结构的单独解释率以及两者的共同作用。以上操作在R软件(4.2.1)中使用vegan包完成。绘图采用软件Origin(2022)、R软件(4.2.1)、Canoco(5)和Adobe Photoshop CC(20.0.7)。

2 结果与分析 2.1 环境因子特征

双因素方差分析显示: 大多数环境因子在季度之间存在显著性差异(如硝态氮、氨氮等)(P＜0.05)(表 1)。枯水期的氮盐、溶氧、矿化度、盐度和氧化还原电位最高, 平水期的总磷、高锰酸盐指数、流速、河宽和水深最高, 丰水期拥有最高的气温、水温、电导率和pH。仅有少数环境因子在河段之间差异显著(P＜0.05)(如总磷、电导率等)(表 2)。上游具有最高的海拔、氧化还原电位和亚硝态氮, 中游具有最高的硝态氮、电导率、矿化度、溶氧、pH和盐度; 下游具有最高的总氮、总磷、氨氮、气温、水温、河宽、流速和水深。此外, 季度和河段的交互作用仅在高锰酸盐指数和气温中差异显著。

2.2 福寿螺种群结构时空特征 2.2.1 密度时空特征

在51个样点中, 有32个样点采集到福寿螺, 出现率为62.7%。共采集福寿螺1771只, 平均密度(2.21±2.35) 个/m²。双因素方差分析显示: 福寿螺密度与季度(枯平丰)、河段(上中下游)以及它们的交互作用均无显著差异(P＞0.05)。在时间尺度上呈现出丰水期((2.57±2.63) 个/m²)>平水期((2.12±2.25) 个/m²)>枯水期((1.88±2.02) 个/m²)的趋势(图 3), 在空间尺度上呈现中游((2.85±2.81) 个/m²)>下游((1.83±2.19) 个/m²)>上游((1.66±1.16) 个/m²)的趋势。

2.2.2 可量性状时空特征

双因素方差分析显示: 福寿螺的体重、壳高、壳宽、壳口宽和壳口长在季度(枯平丰)、河段(上中下游)以及它们的交互作用均在显著差异(P＜0.01)。在时间尺度, 丰水期福寿螺具有最大的体重((12.42±10.15) g)、壳高((35.66±11.40) mm)、壳宽((31.15±10.41) mm)、壳口宽((17.88±6.17) mm)和壳口长((25.75±8.35) mm), 而枯水期各性状具有最小值((9.27±9.33) g、(30.29±12.08) mm、(26.43±11.04) mm、(15.06±6.24) mm、(21.77±8.85) mm)。在空间尺度, 上游福寿螺具有最大的体重((13.32±9.87) g)、壳高((35.81±12.12) mm)、壳宽((31.26±11.13) mm)、壳口宽((17.94±6.57) mm)和壳口长((25.00±8.98) mm), 而下游各性状具有最小值((8.50±8.61) g、(29.77±11.48) mm、(25.99±10.55) mm、(14.87±6.24) mm、(21.56±8.45) mm)(表 3)。

2.2.3 种群结构时空特征

三个季度共采集幼螺694只, 中螺960只, 成螺117只, 年龄结构强烈向中螺倾斜。其中A—J分别占总数的16.1%、11.1%、12.0%、14.5%、14.6%、14.3%、10.8%、4.5%、1.6%、0.5%(图 4)。

总体而言, 福寿螺中螺占比最高(54.2%), 年龄结构为稳定型。双因素相似性分析(two-way ANOSIM)表明: 福寿螺种群结构在时间尺度上(R=0.025, P=0.21)和空间尺度上(R=0.014, P=0.26)的差异均不显著。在时间尺度, 幼螺占比最高, 为49.6%, 年龄结构为增长型; 平水期和丰水期中螺占比最高, 分别为53.2%和61.5%, 年龄结构为稳定型。在空间尺度, 中螺在上游和中游占比最高分别为61.5%和55.2%, 年龄结构为稳定型; 下游幼螺占比最高为48.9%, 年龄结构为增长型(图 5)。

2.2.4 性别比例时空特征

三个季度共采集雌螺1426个, 雄螺345个, 平均性比为0.36±0.40, 其中雌性比例为80.5%, 雄性比例为19.5%。双因素方差分析显示: 福寿螺性比在季度、河段及它们的交互作用均无显著差异(P＞0.05)(表 4)。性比在不同季度呈现丰水期(45.0%)>枯水期(29.1%)>平水期(25.3%)的趋势, 在不同河段呈现中游(35.8%)>下游(34.4%)>上游(28.1%)的趋势。

2.3 福寿螺种群结构驱动机制

CCA结果表明福寿螺种群结构受到环境和空间因子的共同作用, 但不同季度的驱动因子存在一定差异(图 6)。影响枯水期福寿螺种群结构的关键因子为pH、海拔、溶氧、亚硝态氮、盐度、PCNM9和PCNM12。CCA第一轴与第二轴累积解释54.4%的种群结构变异。其中幼螺与海拔、亚硝态氮和溶氧成正相关, 与盐度、pH和PCNM12成负相关; 中螺与溶氧、盐度、PCNM9和PCNM12成正相关, 与海拔和pH成负相关; 而成螺与盐度、pH和PCNM12成正相关, 与海拔、溶氧、亚硝态氮和PCNM9成负相关。

影响平水期福寿螺种群结构的关键因子为流速、矿化度和PCNM5。CCA第一轴和第二轴累积解释41.1%的种群结构变异。其中幼螺与流速和PCNM5成正相关, 与矿化度成负相关; 而中螺、成螺与矿化度成正相关, 与流速和PCNM5成负相关。

影响丰水期福寿螺种群结构的关键因子为水深、高锰酸盐指数、PCNM2、PCNM8、PCNM9和PCNM11。CCA第一轴和第二轴累积解释了52.6%的种群结构变异。其中幼螺与PCNM2、PCNM8、PCNM9和PCNM11成正相关, 与水深和高锰酸盐指数成负相关; 而中螺、成螺与水深和高锰酸盐指数成正相关, 与PCNM8、PCNM9和PCNM11成负相关。

变差分解显示, 福寿螺种群结构受到环境和空间因子的共同影响, 但不同季度驱动因子的解释率存在一定差异(图 7)。三个季度总解释率介于13%—24%之间。此外, 除枯水期外, 其它季节的环境因子解释率均大于空间因子解释率。

3 讨论 3.1 福寿螺种群结构时空格局

福寿螺是典型的外来入侵种。有关福寿螺的研究多见于中小型水体(如池塘、沟渠等)^[21], 而流域尺度的调查极为不足。桂江流域福寿螺密度为(2.21±2.35) 个/m², 该密度明显低于段清星等^[9]在静水池塘中的研究结果((18.14±3.87) 个/m²)。栖息地类型是影响福寿螺密度的重要因素^[22]。不同栖息地包含着不同的食物资源和生境条件等, 中小型水体的营养物质更为聚集、环境稳态更高, 更利于福寿螺的生长和繁殖^[23], 而桂江流域福寿螺受限于多种外界环境压力(如过快的流速导致营养物质贫乏, 且不利于福寿螺附着), 从而具有较低的密度^[24]。其次, 本研究发现福寿螺密度在季度间无显著差异, 表明该物种具有较宽的生态位, 可适应不同季度的环境变化。然而, 福寿螺密度仍呈现出丰水期>平水期>枯水期的趋势, 这与管强^[25]的研究结果一致。丰水期温度高、光照周期长, 附着藻类大量生长, 丰富的食物资源和适宜的温度促进了福寿螺的摄食, 间接影响其产卵行为^{[12, 26]}。此外, 福寿螺密度在河段间无显著差异, 表明福寿螺已扩散至桂江全流域。黄红英等^[23]发现福寿螺在人类影响较大的河流分布较广, 而本研究发现福寿螺在干扰较小河流也有分布, 亟需引起人们重视。因此, 本研究建议不仅在福寿螺广泛分布的季度和河段采取防治手段, 在分布较少的地区也应加强对福寿螺的监测和防治工作, 根据福寿螺种群时空分布规律, 从薄弱处着手, 及时控制灾情。

福寿螺所在栖息环境决定其获得的资源和受到外部环境的干扰程度, 因此栖息环境密切地影响福寿螺的可量性状^[27]。本研究发现福寿螺可量性状在丰水期和上游具有最大值, 而在枯水期和下游值最小, 这与房苗等^[28]研究结果相似。福寿螺在不同干扰程度下的能量分配方式不同。例如, 在干扰程度较低的上游, 福寿螺将更多能量用于自身生长, 进而增强其对环境耐受性, 并扩大对新生境的入侵性和适应性, 表现出较低的密度和较大的个体, 而中下游福寿螺受到外部环境(如天敌、人为干扰等)影响较大, 会倾向于将能量投入到繁殖中, 以延续种群发展, 从而表现出较高的密度和较小的个体^[29]。由此可见, 福寿螺在不同干扰程度下的性状差异和能量分配方式均体现出福寿螺应对外界胁迫的生长—繁殖权衡机制。丰水期温度((30.23±2.43)℃)处于福寿螺最佳温度范围内(25—32℃), 更适宜福寿螺生长发育, 因而有最大的个体, 而平水期和枯水期气温较低, 抑制福寿螺的生长发育, 因此很难长成较大的个体^[12]。

年龄结构可以直接影响出生率和死亡率, 因此种群数量的变化趋势可以通过种群的年龄结构来预测^[30]。本研究发现福寿螺年龄结构整体为稳定型, 但在不同时空中存在差异(均处于增长型或稳定型), 这与房苗等^[28]研究结果类似。桂江流域处于亚热带季风区, 全年平均气温较高且食物丰富, 为福寿螺生存提供了良好条件^[6]。其次, 福寿螺繁殖力明显高于其它螺类, 1年可发生2—3代, 幼螺3—4个月即达到性成熟, 且世代重叠^[31]。因此适宜的栖息环境和强的繁殖力是导致福寿螺种群结构呈现增长型或者稳定型的主要原因。

性别比例在一定程度上影响种群的出生率。本研究发现桂江流域福寿螺性比在季度与河段间均无显著性差异, 且性比通常小于1, 这与Yusa Y等^[32]研究一致。雌螺在资源和空间的竞争中占据了较大的优势, 雌螺生长速度明显快于雄螺, 这得益于雌螺有着更高的摄食率、更强的食物转化率和耐受能力^[33]。同时一些研究表明雌性比例的升高在一定程度上促进种群的数量增长, 雌多雄少有利于提高生殖力和延续种群发展^[34]。

3.2 福寿螺种群结构驱动机制 3.2.1 种群结构与环境、空间因子的关系

本研究初步甄别了显著影响福寿螺种群结构的环境因子(如pH、盐度等)和空间因子(如PCNM2、PCNM5等), 表明环境和空间因子共同影响福寿螺种群结构。

pH、盐度、海拔和溶氧是影响枯水期福寿螺种群结构的主要因素。朱丽霞等^[26]发现pH过酸或过碱均会影响福寿螺的存活, 并限制其体重增长, 降低其产卵的数量及卵的孵化率。本研究区域的pH处于福寿螺的生态范围内, 随着pH的升高, 成螺数量增加, 表明弱碱性的水体更利于成螺的存活和碳酸钙外壳的形成, 这与Seuffert等^[35]研究结论一致。海拔和盐度呈现出相反的变化趋势, 两者均是影响福寿螺分布的重要因素^[36—37]。本研究发现幼螺多分布在海拔较高、盐度较低的上游, 表明幼螺的生长发育在高海拔、低温的环境下受到抑制, 同时幼螺活动能力差, 易于捕捉^[38—39]。而海拔较低的中下游地区, 温度适宜, 合适的盐度又能为福寿螺提供更加丰富的无机盐离子^[36], 因此, 该区域主要以中、成螺为主。水体溶氧亦是影响水生生物分布的重要因素, 低氧会明显影响螺类的生长及分布, 厌氧则会导致其大量死亡^[40]。本研究发现幼螺和中螺密度随着水体中溶氧的上升而增加, 这与高翼天等^[41]结果类似, 表明充足的氧气是促进幼螺生长发育和福寿螺孳生分布的必要条件之一, 而成螺的呼吸管道更为发达, 其耐低氧的能力更强。

流速和矿化度是影响平水期福寿螺种群结构的主要因素。福寿螺分布对流速亦有一定要求。本研究发现幼螺在流速较快的环境中分布较广, 较快的流速不仅保证了充足的氧气供给和丰富的营养物质, 同时也加快了水体更新, 保证了水质的清洁, 有利于幼螺生长发育^[42], 这与高翼天等^[41]研究结果类似。矿化度可通过影响水生植物的群落组成从而间接影响福寿螺种群结构^[43]。本研究发现中螺和成螺在矿化度较高的环境中分布较广, 其原因可能是较高的矿化度有利于水生植物(如苦草和黑藻)的生长, 丰富的水生植物为福寿螺提供了充足的食物资源; 福寿螺为广食性种类, 中螺和成螺主食鲜嫩水生植物, 其种群结构分布与摄食的植物丰富度有关, 这与龙伦明等^[43]结果类似。

高锰酸盐指数和水深是影响丰水期福寿螺种群结构的主要因素。高锰酸盐指数反映水体中还原物质含量, 代表水体污染程度。有研究报道, 福寿螺偏爱水草丰盛、水质清新的水域, 但也能忍受污浊的水质^[44]。本研究发现中螺和成螺在水质较差的环境分布较广, 这可能是由于中螺和成螺抗逆性功能较为完善, 而发育初期的幼螺抗逆性低。水深是福寿螺分布的重要限制因子之一, 水深的变化会引起一些重要指标(如溶解氧和水温等)的改变, 从而间接影响福寿螺的种群结构^[45]。本研究发现幼螺常栖息在较浅的水体中, 其原因可能是较浅水域中有大量泥土基质, 当遇到不良环境时, 其能迅速躲避, 同时浅水区的食物资源较丰富, 不仅为幼螺提供了生长发育的营养物质, 还更有利于其进行水生和陆生的交替生活, 这与前人研究结论一致^[46]。

此外, 空间因子所表征的扩散过程同样影响着福寿螺的种群结构, 且在不同时期存在一定差异, 但均为中、大尺度的空间变量(枯水期: 中尺度; 平水期: 大尺度; 丰水期: 中、大尺度)。不同尺度的空间因子决定福寿螺是否有机会到达最适的生境进行定殖、生长和繁殖, 最终改变福寿螺的种群结构。扩散过程的重要程度与不同年龄种群的扩散能力、各生境间的距离和连通度有关^[47]。

3.2.2 环境、空间因子对种群结构的相对作用

本研究发现不同季度驱动力的相对作用存在差异, 这可能与环境因子的季节性变化以及福寿螺生活史有关。空间因子对枯水期福寿螺种群结构的作用更为显著。福寿螺营两栖生活, 主要借助水环境进行入侵和扩散^[12]。枯水期温度降低, 福寿螺的活动能力下降, 同时, 桂江流域为树枝状的水系, 河网复杂, 且采样点之间的距离相隔较远, 多种因素阻碍了福寿螺的扩散, 意味着此时空间作用(扩散限制)的影响比较强烈^[48]。因此枯水期是消灭福寿螺、降低繁殖基数及防止来年灾害爆发的最佳时期。

环境过滤对平水期和丰水期福寿螺种群结构作用更为显著, 其强度受到环境因子的限制, 导致不同年龄结构分布在特定的环境梯度^[49]。由于河流生态系统的复杂性和异质性, pH、海拔和溶氧等理化特征存在着一定的变化梯度, 同时平水期和丰水期水位上升, 河流连通性增强, 福寿螺可沿着河道扩散寻找最适的生境进行定殖。环境因素作为过滤器, 使得河流生态系统中环境过滤的作用往往占上风^[50]。因此应针对不同发生区域, 因地制宜、分类施策, 开展防治。加强未发生区域的检疫监控和防控力度, 在螺害发生较重的地区, 设置拦截网捕螺灭螺, 避免其扩散繁殖。

本研究中采用的变差分解分析方法并不能解释福寿螺种群的所有变异, 表明其它因子同样起着重要的作用, 其中包括人类行为导致的环境破碎化(如水生植物的消亡和河岸带固化)^[51]、生物间相互作用(如种间竞争、天敌捕食)等^[24]。

4 结论

本文初步探究了桂江流域福寿螺种群结构分布特征及驱动机制。结果表明桂江流域福寿螺种群现已广泛分布整个流域, 并且种群处于稳定期(以中螺为主)。福寿螺的密度在丰水期和中游具有最大值, 在枯水期和上游值最小; 而可量性状在在丰水期和上游具有最大值, 在枯水期和下游值最小, 表明该螺在应对不同环境条件时能量分配方式可能有所不同。同时, 桂江流域福寿螺雌多雄少的性别比例表明该种群依然存在进一步扩散的风险。最后, 福寿螺种群结构驱动机制在不同季度间存在差异, 枯水期主要是受扩散限制的影响, 而平水期和丰水期主要受环境过滤的影响。