水库具有特殊的水文和生态学特征, 其生态系统状态的变化可能威胁到生物多样性, 同时对居民用水安全造成重大影响[1, 2].在水库运行过程中, 营养物质持续地在流速相对较缓的入库水体的携带下进入水库, 并在库区中滞留、沉降及富集, 最终导致水体富营养化[3].近年来, 水库面临的富营养化风险日益严重, 如三峡水库[4]、广西大龙洞水库[5]和湖北陆水水库[6]、俄罗斯Ivankovo水库[7]、美国Eagle Creek水库[8]和巴西Funil水库[9]等均表现为富营养化状态.城市用水主要来源为水库, 如果不对水库富营养化过程加以干预、控制, 将给自来水厂及居民用水安全等造成难以估量的危害.因此, 密切监测水库水质变化, 对其生态系统演变特征进行探究迫在眉睫.

在对水库系统中的营养物质处理机制研究中, 浮游微生物常被认为是元素及营养物质转移的主要参与者, 同时被视为能量流动的重要贡献者[10~12].在以往的研究中, 研究者主要对水库中的细菌群落进行了广泛且深入地探究[2, 5, 9, 13, 14].但是浮游微生物在生态功能中的作用不仅归因于细菌群落, 还归因于微生物生态系统的其他组成部分, 如真核微生物群落[1, 15].有研究表明, 真核微生物是水生食物网的关键组成部分, 具有重要的生态功能, 其既是有机质的生产者, 也是细菌生物量的主要消费者[16, 17].淡水生态系统中真核微生物群落的动态变化与环境变量密切相关[18], 并且有研究发现富营养化水域的物种多样性较少[19].有研究发现, 真核微生物分布广泛且种类繁多, 对维持生态平衡和促进水生生态系统中的物质及能量循环至关重要[15].但是与细菌等原核微生物相比, 饮用水水库中真核微生物的群落组成和动态研究相对较少[1].对水库真核微生物群落组成进行研究, 了解其在水库中的分布特征及其在水库生态系统中的作用, 对保护和维持水库生态系统的平衡具有重要意义[20].

目前, 对水体中的真核微生物的研究主要通过显微镜观察, 但是这种方法不仅受到采样条件和保存技术的限制, 而且在鉴定物种类别时也存在很大的误差[20].随着高通量测序技术的出现和发展, 可以更为高效快速地获取特定的DNA片段, 并且得到更全面地群落组成[21], 更准确地评估真核微生物多样性.因此, 高通量测序被广泛地应用于水体中真核微生物群落结构, 进化关系以及真核微生物与环境之间相关性的研究中[20, 22], 以探究真核微生物的群落组成, 分布特征和影响因素[23, 24]. Keck等[25]的研究利用高通量测序技术对48个湖泊中微真核生物群落进行对比发现, 其群落组成存在显著的差异, 并且在全球范围内, 偏远的高海拔湖泊受当地人类活动的影响小于低海拔湖泊, 同时研究表明人类活动的变化有利于光营养型和混合营养型物种的生态位竞争. Reich等[24]的研究发现, 原位海水和pH降低调整后的海水中真菌群落结构和多样性存在显著差异, 主要表现为优势群落间的丰度变化; 同时在pH值较低的海水中, 微生物群落中占主导地位的一些真菌是已知的病原真菌.共生网络分析方法是一种强大的生物信息分析技术, 已广泛应用于探索水库[26]、湖泊[27]和海洋环境[28]中微生物群落间复杂的相互作用关系.Yan等[1]的研究利用共生网络分析发现金盆水库中真核生物对浮游植物群落动态结构的变化具有一定响应, 同时两者之间可能存在拮抗或互惠关系. Zheng等[20]的研究通过共生网络分析发现丹江口水库真核生物群落中的关键类群及影响其结构的关键水质因子.但是, 关于水源水库真核微生物群落季相演替特征和驱动因素的研究却鲜见报道.

鉴于此, 本文采用高通量测序技术, 以陕西省西安市李家河水库为研究对象, 在对水库水质跟踪监测基础上, 通过调查研究李家河水库真核微生物群落结构组成和水质理化条件变化, 结合共生网络分析对不同真核微生物种群间的共生与相互作用, 及其与水质因子间的响应关系进行解析, 明确了真核微生物的群落结构季相演替特征和驱动因素.具体研究如下: ①探明李家河水库中的主要水质参数的变化; ②通过高通量测序技术研究真核微生物群落结构的季相演替规律; ③通过共生网络分析和多元数据统计方法解析真核微生物不同种群间的复杂联系及影响其群落结构的驱动因素, 以期为研究水生生态系统中的真核微生物和环境变量提供科学依据.

李家河水库位于陕西省西安市蓝田县南面辋川镇境内, 距离县城大约19 km, 库区位于33°58′30″~33°59′41″N, 109°23′36″~109°24′36″E之间.该水库于2015年运行, 是一座以向西安市浐河以东区域供水为主, 兼具发电和防洪功能的水库, 总库容5 690万m3 [13].随李家河水库的持续运行, 该水库面临着富营养化的潜在风险[29].

本研究在李家河水库设置S1(33°59′28″N; 109°23′46″E)、S2(33°59′17″N; 109°24′00″E)和S3(33°59′34″N; 109°23′55″E), 共3个采样点.采样工作进行于2020年5月至2020年10月期间, 每月于采样日下午2时~6时在各采样点2 m深度处采集水样2.5 L, 用于部分水质参数总氮(TN)、总磷(TP)、氨氮(NH4+-N)、硝氮(NO3--N)、总溶解性有机碳(DOC)、铁(Fe)、锰(Mn)和叶绿素a(Chl-a)的测定以及真核微生物群落结构高通量测序.温度(T)、pH、溶解氧(DO)、电导率(conductivity, Cond)和浊度(turbidity, Tur)使用Hydro-lab DS5型(HACH, USA)多参数水质分析仪进行原位测定[2].水质参数均参照文献[30]进行测定, 所有水质参数分析在48 h内完成(n=3).

水样(1 L)用已灭菌的0.22 μm聚碳酸酯膜(Millipore, 美国)过滤, 用于真核微生物的富集[1].将滤后膜粉碎并转移至已灭菌的50 mL离心管中, 置于低温冰箱(-20℃)保存.样品中的总DNA通过DNA提取试剂盒(Water DNA Kit, Omega, 美国)提取[20], 而后使用Nano Drop 2000(Thermo Scientific, 美国)确定其浓度和纯度[14].

利用真核微生物特异性引物1380F_1510R扩增真核微生物18S rRNA基因的高变区V9区, 从不同样品中获得目标DNA序列[31], 每个样本进行3次重复.PCR扩增[20]采用TRAGEN AP221-02, TransStart Fastpfu DNA聚合酶体系, 混合液(20 μL)由5×反应缓冲液(4 μL)、2.5 mmol·L-1 dNTPs(2 μL)、5 μmol·L-1正向引物和反向引物(0.8 μL)、快速Pfu DNA聚合酶(0.4 μL)、BSA(0.2 μL)、模板DNA(10 ng)和ddH2O组成; PCR反应包括95℃下初始变性3 min, 95℃下变性30 s, 55℃退火30 s, 72℃延伸45 s, 72℃延伸10 min, 共35个循环.通过Illumina MiSeq PE300(Illumina, 美国)测序平台对真核微生物群落结构组成进行分析[26], 测序工作在上海凌恩生物有限公司进行.

本研究中所有水样中真核微生物18S rRNA基因测序结果均已上传至NCBI-SRA(美国国立生物技术信息中心)数据库中(网址: http://www.ncbi.nlm.nih.gov/), 编号为SRP297316和SRP309745.

使用Excel 2019对本研究中所有数据进行初步梳理.水质参数时间变化单因素方差分析(One-way ANOVA)使用软件SPSS(version 24.0)进行处理, 显著性水平的参数选取0.05, 并进行Tukey HSD事后检验.在QIIME 2(微生物生态学定量分析, version 1.9.1)中将小于50 bp的真核微生物序列读数进行过滤预处理, 然后将相似性较高(≥97%)的序列归为一个OTU分类单元[26].使用MOTHUR软件进行真核微生物种群多样性计算.分别使用Circos网站在线绘图功能(http://circos.ca/)和R语言绘制热图(Heatmap)对研究期间真核微生物种群门水平和属水平上的相对丰度变化情况进行可视化表达.为了更直观深入地探明不同真核微生物之间的互作规律以及其同各水质指标之间的变化响应关系, 本研究将使用Gephi软件(version 0.9.2)对藻类和各水质指标绘制共生网络图, 参数设置为P < 0.05, |R|>0.8[28].此外利用CANOCO 4.5软件对真核微生物与水质因子的关系进行RDA分析.

本研究期间李家河水库各水质参数随时间动态变化情况如表 1所示, 从春末5月至夏初6月水温大幅上升, 而后在整个夏季(6月~8月)及秋初9月保持相对较高水平, 并在8月达到最高水温(20.35±0.64)℃, 而后进入秋季开始呈下降趋势.DO随着水温的上升呈现相反的变化趋势, 在夏季含量较低, 而在水温较低的春末及秋初的水体中具有较高的含量.Diaz等[32]的研究表明, 溶解氧浓度可以反映水质的污染程度, 特别是与有机物含量有关的污染程度, 是衡量水质的重要指标.ρ(NH4+-N)在夏初达到最高值(0.31±0.01)mg·L-1, 在7月降到(0.12±0.02)mg·L-1, 随后在夏末及秋季保持相对稳定.ρ(TN)和ρ(NO3--N)均在春末5月达到最高值, 分别为(3.09±0.17)mg·L-1和(2.58±0.31)mg·L-1, 进入夏初两者均有较为明显的下降, 而后在7月上升; 同一时期Chl-a表现出相反变化趋势, 可能是因为藻类的增殖与消亡对水体中氮的利用与释放有关. ρ(TP)在整个研究期间变化不大, 最高出现在秋初9月(0.11±0.02)mg·L-1. ρ(Fe)和ρ(Mn)均小于集中式生活饮用水地表水源地标准限值, 其中ρ(Fe)整体呈现先上升后下降趋势, 在7月达到最高值(0.16±0.01)mg·L-1, 有研究表明较高的水中铁浓度可促进藻类水华的增殖[33], 而藻类与真核微生物之间存在着密切的相互作用[1], 而ρ(Mn)在研究期间变化不大.浊度在整个研究期间处于较低水平, 7月与10月的高浊度值主要是因为在7月与10月采样期间李家河水库受到日径流输入的影响, 故其剧烈变化与径流的输入有密切的关联[34].总体来说, 根据“地表水环境质量标准”(GB 3838-2002), 本研究期间李家河水库的pH符合地表水环境质量标准; DO值均符合地表水年度Ⅰ类标准; ρ(NH4+-N)符合地表水Ⅱ类标准, 而ρ(TP)和ρ(TP)偏高.本研究期间监测结果显示, 李家河水库水质整体较为良好, 与以往研究的结果一致[2, 13].

在本研究中, 通过采用α多样性分析, 对李家河水库不同时期的真核微生物群落结构组成动态变化进行表征, 包括物种丰富度指数(Chao1), Shannon多样性指数(H′)和Simpson多样性指数(D).对李家河水库3个采样点样品测序后, 共获97%相似水平及以上的OTUs有3 230个, 且对真核微生物群落覆盖率达到99%(表 2).OTU的最大值612出现于夏季8月, 而在夏初6月时, OTU数达到全年最小值464.同时可以发现, 本研究中的OTU值于夏季整体较春秋季节高.对研究期间3个采样点样品测序结果进行分析, 表 2显示, 3个采样点的Chao1、Shannon和Simpson指数均值均在研究周期内差异性表现较为明显.其中, Chao1丰富度指数在春季整体表现最低, 其次为秋季, 表明李家河水库在夏季的真核微生物群落丰富度高于春秋季节, 秋季真核微生物群落丰富度又高于春季.此外, 还可以从表 2中发现研究期间每月对应的真核微生物群落OTU值均低于Chao1丰富度指数, 表明在李家河水库中存在许多未知真核微生物序列[18].Shannon和Simpson指数的最大值分别出现于8月(4.14)和10月(0.349), 而两者最小值分别出现于10月(2.35)和6月(0.056), 表明李家河水库真核微生物群落多样性在6月和8月达到最大, 而在10月时多样性最小.综上所述, 李家河水库真核微生物群落丰富度和多样性均表现出夏季高于春秋季.

李家河水库研究期间, 对真核微生物群落门水平选取OTU相对丰度的前19具有代表性的物种进行Circos分析.如图 1所示, 真核微生物群落相对丰度在门水平上有明显季相的变化.从中可以发现, 真核生物物种丰度在夏季时整体高于春秋季节.在春末5月相对优势真核微生物为Kathablepharidae和节肢动物门(Arthropoda), 两者之和达到了54.30%, 其次为褐藻门(Ochrophyta)和Rigifilida分别占比10.81%和9.49%.夏初6月, 节肢动物门(Arthropoda)继续保持其相对丰度优势(23.84%), 同时为优势物种的是纤毛门(Ciliophora, 25.21%)和Phragmoplastophyta(17.65%).有研究表明, 纤毛虫是细菌捕食者之一[19], 同时与藻类间存在共生关系[35], 说明纤毛虫不仅与环境因素有关, 还受到生物因素的影响.到达7月后, 节肢动物门(Arthropoda)成为该月绝对优势物种, 占比达到了60.81%, 而后进入8月衰减至2.83%.王昱等[36]的研究发现节肢动物为底栖动物群落中的重要类群, 并在水质为清洁~轻污染情况下为优势类群.从夏末8月至秋季, 李家河水库水体中的真核微生物种群结构发生了改变, 其中较为显著的是节肢动物门(Arthropoda)不再是优势门类, 而刺胞动物门(Cnidaria)、Others和绿藻门(Chlorophyta)成为新的优势门类, 其中Others在8至10月占比逐月下降, 刺胞动物门(Cnidaria)在秋季10月成为绝对优势门类(77.27%), 绿藻门(Chlorophyta)则是在9月成为优势门类(45.79%).研究表明, 刺胞动物门(Cnidaria)组成和丰度在时间上的变化是常见的, 可能是由局部尺度上的生物因素引起, 例如藻类的结构复杂性[37], 藻类衰亡以及与附生藻类或其他种类的水生动物对底物的竞争[38, 39].

采用热图(Heatmap)对李家河水库研究期间真核微生物属水平的群落结构变化(选取前29具有代表性的丰度物种)进行表征, 如图 2所示.随着时间的推移, 李家河水库真核生物群落组成在属水平上表现出明显的多样性和差异性.在春季末5月, Telonema sp.出现并在当月表现为绝对优势, 其相对丰度为50.53%.Telonema sp.在海洋环境中被广泛地报道[40], Lefèvre等[41]的研究通过高通量测序首次在淡水环境中发现它的存在.进入夏季后, 初期6月扁藻属(Tetraselmis sp.)和隐核虫属(Cryptocaryon sp.)为最丰富的属类, 占比分别为30.11%和24.23%, 另一具有较高相对丰度的属类为Phascolodon sp., 其占比为10.15%, 其余属类相对丰度占比均低于10%.扁藻属多见于海洋环境, 但是也有一些淡水种, Hu等[4]对三峡地区真核浮游植物时空演替规律的研究发现, 扁藻属在夏季为优势属类, 并且有潜在的藻华暴发风险.隐核虫属是一种常见的鱼类寄生虫[42], 其在夏季的大量繁殖可能与鱼类产卵繁殖有关.鱼鳞藻属(Mallomonas sp., 20%)、隐藻属(Cryptomonas sp., 15.79%)、杯鞭虫属(Bicosoeca sp., 15.79%)、Poteriospumella sp.(10.66%)和隐核虫属(Cryptocaryon sp., 10.26%)在夏季中期7月成为相对丰度较高的属类.Maznah等[43]的研究表明, 鱼鳞藻属为水体水质朝不良状态发展的潜在指标.隐藻属广泛存在于不同水体环境中, 并且其没有明显的季节分布特异性[4].但是由隐藻属大量繁殖引发的水华现象也有报道[44].此外, 有研究报道衣藻属和隐藻属在藻类群落中发挥着主导作用[44, 45].夏季末8月真核微生物优势属类为Junceella sp., 占比高达75.51%.在秋初9月, 扁藻属(Tetraselmis sp.), Junceella sp.和衣藻属(Chlamydomonas sp.)为优势真核微生物属类, 占比分别为41.34%, 20.55%和20.31%.10月Junceella sp.又重新成为绝对优势真核微生物属类, 占比高达87.18%.目前, 对于Junceella sp.在淡水环境中的生态学特性鲜有报道.

为探究真核微生物间的互作关系, 本研究采用共生网络分析技术对研究期间李家河水库相对丰度位于前40位的具有代表性的真核微生物属水平的物种进行物种间互作关系分析(图 3).该网络由38个节点及312条边组成. 不同颜色的节点代表不同的类群, 节点大小代表分类学丰度, 绿线代表负相关, 红线代表正相关, 线粗代表相关性的强度, 粗线表示相关性较强[20].平均度大小与节点之间联系紧密性呈正相关, 本网络图 4中平均度为16.421, 表明李家河水库真核微生物之间的联系较为紧密.所有的节点被划分为3个模块, 分别包含13、17和8个节点(属), 模块是指一组彼此连接良好但与其他模块的OTU联系较少的OTU, 在网络中, 同一模块中的大多数OTU(节点)具有密切的相关关系.从图 3中可以发现, 每个模块具有不同的组成特性, 在模块Ⅰ中真核微生物间主要以正相关关系为主, 而在模块Ⅱ与模块Ⅲ中则大多数为负相关关系.同时也可以发现, 不同模块间也存在紧密的联系, 并以负相关关系为主.正相关关系表明不同类群的真核微生物的生态位较为相似或彼此之间存在互利共生的关系, 而负相关则表明他们之间可能生态位不同或彼此之间存在竞争、拮抗的作用[46].Phascolodon sp.仅与斜片藻属(Plagioselmis sp.)间出现负相关性(P < 0.01), 说明两者间存在一定的捕食或者竞争关系, 与这两者类似的有亚历山大藻属(Alexandrium sp.)和弹跳虫属(Halteria sp.).Spumella sp.也仅与扁藻属(Tetraselmis sp.)存在负相关性(P < 0.01), 不同的是其两者均为真核藻类, 两者之间可能存在着营养底物的竞争.杯鞭虫属(Bicosoeca sp.)、尖毛虫属(Oxytricha sp.)和壳虫属(Tintinnidium sp.)则表现出与其他真核微生物间强烈的相关性且以负相关性为主(P < 0.01), 表明其对李家河水库真核微生物群落结构起着关键作用, 且在物种间多存在着竞争关系.最近, Zheng等[20]通过共现网络分析对丹江口水库真核生物群落结构特征的研究发现, 群落结构随时空变化, 并且不同类群间的连通性较高. Xue等[26]的研究使用共现网络分析发现, 在亚热带水库中稀有类群和非稀有类群的协同作用可能在维持真核微生物群落的稳定性和生态功能方面发挥重要作用.

对真核微生物与水质进行网络分析, 网络图 4中共有46个节点, 370条边.由图 4可知, 水质参数中, pH、Chl-a和电导率对真核微生物的影响较大, 且pH和Chl-a对真核微生物的影响以正相关为主(P < 0.01), 而电导率主要呈现出与真核微生物的负相关性(P < 0.01).Chl-a与歪水蚤属(Tortanus sp.)、隐核虫属(Cryptocaryon sp.)和挠足目温剑水蚤属(Thermocyclops sp.)存在负相关关系(P < 0.01), 表明这些真核动物可能与真核藻类间存在一定的竞争关系.同时可以发现, 囊担菌属(Cystobasidium sp.)与水温呈现正相关性(P < 0.01), Sen等[47]的研究发现水温是影响其多样性的重要因素.此外, 有研究表明, DO、Chl-a、水温、pH和电导率可能是影响浮游真核生物的主要因素[48].

对李家河水库真核微生物群落与水质关系进行冗余分析(RDA), 结果表明(图 5), 在研究期间, 李家河水库真核微生物群落结构受到不同水质的影响具有显著差异.水质共解释了真核微生物群落组成的86%的变异, 其中RDA1解释了69.4%的变异, RDA2解释了16.6%的变异.从图 5中可知, 影响真核微生物群落结构最重要的决定因素依次为浊度、pH、DOC、TN、DO、TP、水温和电导率(P < 0.05).Chen等[49]的研究发现, 真核微生物群落结构主要与pH、DO和DOC浓度的变化有关, 这与本研究的结果基本一致.不同时期, 水质对真核微生物群落组成的影响也具有一定的差异, 真核微生物群落结构组成在春季末及夏季主要受到TN、电导率、pH、DO和Mn浓度的影响(P < 0.05), 而在秋季TP、DOC、水温和浊度与其具有更为紧密的联系(P < 0.05).

(1) 对研究期间李家河水库水质进行连续监测结果表明, 水库水质总体良好, pH处于(7.49±0.13)~(9.26±0.14), ρ(DO)最低为(7.80±0.33)mg·L-1, ρ(NH4+-N)最高为(0.31±0.01)mg·L-1, 整体符合Ⅱ类地表水标准; ρ(Fe)最高为(0.16±0.01)mg·L-1, ρ(Mn)最高为(0.03±0.00)mg·L-1, 两者均小于集中式生活饮用水地表水源地标准限值, 但总氮和总磷含量仍较高.

(2) 对高通量测序所获的97%相似水平及以上的3 230个OTUs进行多样性分析结果表明, OTU值于夏季整体较春秋季节高, 多样性指数在研究期间呈现出较为显著的动态变化特征, 真核微生物群落Chao1丰富度指数和Shannon多样性指数(H′)均表现出夏季高于春秋季.

(3) 在门水平分析表明, 李家河水库真核生物群落组成具有一定的季节性差异.节肢动物门在春末至夏季中期为优势门类.此外, 纤毛门(25.21%)为夏初真核微生物优势门类.刺胞动物门和绿藻门为夏末至秋季优势门类, 而刺胞动物门为秋季绝对优势门类(77.27%).

(4) 李家河水库真核生物群落组成在属水平上表现出较门水平更为明显的多样性和时间差异性.Telonema sp.为春末5月优势属类(50.53%).扁藻属(30.11%)为夏初6月优势属类.鱼鳞藻属(20%)在夏季中期(7月)成为相对丰度较高的属类.Junceella sp.在夏末8月及秋季前中期均有较高相对丰度, 最高占比达87.18%.扁藻属(41.34%)为9月优势属类.

(5) 杯鞭虫属、尖毛虫属和壳虫属表现出与其他物种间存在强烈的相关性且以负相关性为主.RDA分析表明, 真核微生物群落结构组成在春季末及夏季主要受到TN、电导率、pH、DO和Mn浓度的影响, 在秋季主要受到TP、DOC、水温和浊度的影响.网络分析表明, 真核微生物多与pH和Chl-a呈现正相关关系, 与电导率呈负相关关系.