成土时间对黄河三角洲植物及微生物演替的影响

侯冠群,来剑斌,李 静,刘 振,公华锐,王 兵,孙志刚,欧阳竹,侯瑞星, *

1 中国科学院地理科学与资源研究所,生态系统网络观测与模拟重点实验室, 北京 100101

2 中国科学院大学资源与环境学院, 北京 100049

3 中国地质调查局烟台海岸带地质调查中心, 烟台 264004

滨海新生土地演替的研究有助于湿地资源的保护,为改善当地的生态环境提供科学依据[1—2]。了解新生土地生境更新、演替规律,还有助于评估人类活动的影响[3],使开发速度低于生境更新、演替的速度,或者引导和控制新生土地生境更新、演替的速度和方向[4]。

演替的监测需要频繁连续的时间序列。但是,初生演替的研究往往受到时空因素干扰或限制,通常采用时空替代法拟合时间序列,如冰川退缩迹地研究中通常采用样线法[5],而土壤年龄的驱动力研究多采用样地法[6]。在滨海三角洲演替的研究中,采样和分析方法通常为“空间换时间”[7—8]或针对植被演替阶段的样地法[9—11]。虽然前人很早关注到河口三角洲土壤年龄对生态系统演替阶段可能带来的影响[8, 11—12],但缺少成土时间的驱动力的定量研究,这可能造成演替时间和空间的不匹配,进而形成对演替过程的描述产生偏差。

黄河三角洲自1855年黄河夺大清河流路入渤海以来,黄河尾闾多次改道,形成多个呈叠瓦状复式结构的亚三角洲堆积体[13]。其中有4个亚三角洲堆积体与其他洲体的重叠程度较低,在地表形成相对独立的冲积扇区,成土时期分别为1904—1929年,1929—1934年,1964—1976年和1976年至今[14—15]。亚三角洲堆积体重叠时,表层土壤成土时期由晚近流路的冲淤决定[16],故相同亚三角洲堆积体的土壤具有相对一致的发育年龄。各扇区均从内陆向滨海延伸,均包含滨海典型植被群落和植被演替过程,并且均在内陆深处形成了农田土壤。

黄河三角洲土壤含盐量较高且以裸地为起点[17]。随着演替进行,植物种类增加,物种多样性上升,植被盖度、物种丰富度增加[3, 8],土壤质量呈现逐渐改善的趋势[7—8, 18]。植被演替和土壤演替在滨海生态系统内交替发生,且均受到成土时间和离海距离等因素影响。成土时间的不同带来生态系统年龄的差异,且与离海距离共同影响土壤脱盐程度和熟化特征[8, 18]。盐度影响海岸和河口的碳源和汇功能[19],进而影响演替过程。白茅是轻度耐盐植物,可作为可垦殖的植被指标[8],耕地主要由白茅草甸和芦苇-荻草甸开垦而来[3]。黄河三角洲土壤和植被随时间的演替过程也伴随着微生物的演替,微生物群落的演替受土壤化学性质和植物种类影响[10, 12, 17],而微生物群落作为生物地球化学循环的关键因素,又塑造了环境特征和功能[20—22]。

黄河三角洲多个亚三角洲堆积体可提供对滨海土壤向农田土壤演替的时空要素开展研究的平台。本文探究成土时间是否为黄河三角洲植物和微生物群落演替的生态驱动力,为提高评估滨海生态系统植物和微生物演替规律的准确性提供科学依据。

1 材料与方法

1.1 研究区概况

本研究采样区域位于山东省东营市黄河三角洲地区(37°31.807′—38°07.664′N,118°15.858′—119°13.384′ E),地跨东营市河口区、垦利区、利津县。黄河三角洲属暖温带季风气候,年平均气温12.3℃,年平均降水量542.3mm。土壤以滨海盐土和滨海潮土为主, 土壤盐渍化严重, 主要是氯化物盐土和氯化物潮化盐土, 土壤组成以泥沙为主, 养分含量低[23]。自然植被以盐生植物为主。植被带垂直于河床分布,沿河流由海向陆发展[9]。

1.2 样品采集

2020年5月,在各时期冲积扇区范围内进行土样和植被信息采集(图1),按成土历史分别命名为路线1扇区、路线2扇区、路线3扇区和路线4扇区(下文分别简称R1、R2、R3和R4)。分别由河口所在海岸沿河道或古河道方向设计采样路线,每条采样路线由沿海到内陆对样地编号为样地1、样地2……样地6(下文分别简称S1、S2……R6)。采样标准为自成土以来受人类活动干扰较小的、有代表性的原生景观植被,包括典型的盐滩、盐生草地、中生草地和灌丛草地,并采集农田样地。基于实地调查,选取样地21处,如表1所示。R1扇区受人类近海工业和养殖活动影响,未发现原生盐滩；R4扇区成土时间较晚,未发育形成灌丛草地。R1S4、R2S4、R3S5、R4S4为各冲积扇区发现的内陆顶级景观植被,农田样地在以上样地5—10km外的内陆选取。

图1 研究区域和采样地点示意图

表1 样地信息

每处样地采用五点取样法,样点间距5—10m。使用带刻度的土钻在每个样点采集20cm深度土样,将5个样点各深度的土样分别于自封袋中混合均匀。在每个样地的中心样点划出1m×1m样方,在自然群落采集所有草本植物地上部茎叶用于生物量测定,并记录样方内物种信息；记录该处自然群落的各植物类群的盖度信息。

各样地的5个样点的土壤样品混合,分成两部分。其中一部分风干并过20目标准筛,再从其中取约20g过100目标准筛,用于土壤理化性质测定。另一部分立即过20目筛,填充5mL离心管并标号,埋藏于干冰中冷冻,每个样地设三个重复,送至上海美吉生物医药科技有限公司测序。

1.3 土壤理化性质测定

取过20目筛的干土,使用去离子水按水土比5∶1制备浸提液,以PHS—3E酸度计测定pH,以Mettler Toledo电导仪(FiveEasy Plus FE20)测定电导率(EC)。土壤含水量(SWC)用铝盒烘干法测定；土壤有机质(SOM)依据NY/T1121.6—2006标准滴定；全氮(TN)使用浓硫酸消解—凯氏定氮法测定；全磷(TP)和全钾(TK)依据DZ/T0219.1—2006标准使用ARL 9900X射线荧光光谱仪测定；碱解氮(AN)依据LY/T 1228—2015标准滴定；速效磷(AP)依据LY/T1232—2015标准使用TU—1810紫外可见光分光光度计滴定；速效钾(AP)依据NY/T890—2004标准使用ICE3500原子吸收分光光度计测定。以上10个指标连同碳氮比(CN.ratio)作为本研究的土壤理化性质因子集合。

1.4 植被多样性和时空因子计算方法

植被样品取回后,于烘箱中105℃处理24h,称得生物量干重,记为plant.bioM。

植物物种丰富度为样地中所能识别的物种种类数,记为plant.S。统计植物物种数时,考虑到5月许多植物处于苗期或营养生长期而不易辨别,故将同一科的相同生活型的部分物种进行了合并。样地区域植被总盖度记为plant.cov。

植物物种多样性计算方法为Shannon-Wiener指数[24],计算公式为：

H=-∑PilnPi

式中,Pi指第i个物种的出现频率。下文简称Shannon指数。

以样方内植被各物种的多度与植株个体总数的比值为Pi计算物种组成Shannon多样性,记为Shannon.abun。考虑到作为环境因子时,1m2的样方对样地植被信息代表性可能不足,故另以样地植被各类群的盖度占比为Pi计算类群盖度Shannon多样性,记为Shannon.cov。以上5个指标作为本研究的植物多样性因子集合。

各冲积扇区成土时间计算方式为亚三角洲体形成年份区间的中值与采样时间(2020年)的年代差；河口海岸距离计算方法为样地与黄河河口或古河口所在海岸方向的距离。

1.5 微生物群落组成测序方法

根据 E.Z.N.A. soil DNA kit (Omega Bio-tek, Norcross, GA, U.S.)说明书进行微生物群落总 DNA 抽提。使用1%的琼脂糖凝胶电泳检测DNA的提取质量,使用NanoDrop2000测定DNA 浓度和纯度；使用338F (5′-ACTCCTACGGGAGGCAGCAG-3′) 和806R (5′-GGACTACHVGGGTWTCTAAT-3′) 对 16S rRNA基因V3-V4可变区进行PCR扩增,使用SSU0817F (5′-TTAGCATGGAATAATRRAATAGGA-3′)和SSU1196R (5′-TCTGGACCTGGTGAGTTTCC-3′)对18S rRNA基因V4—V5可变区进行PCR扩增。随后基于Illumina MiSeq进行16S rRNA、18S rRNA扩增子测序。使用Fastp(https://github.com/OpenGene/fastp,version 0.20.0)软件对原始测序序列进行质控,使用FLASH(http://www.cbcb.umd.edu/software/flash,version 1.2.7)软件进行拼接：(1)过滤reads尾部质量值20以下的碱基,设置50bp的窗口,如果窗口内的平均质量值低于20,从窗口开始截去后端碱基,过滤质控后50bp以下的reads,去除含N碱基的reads；(2)根据PE reads之间的重叠关系,将成对reads拼接成一条序列,最小重叠长度为10bp；(3)拼接序列的overlap区允许的最大错配比率为0.2,筛选不符合序列；(4)根据序列首尾两端的条形码和引物区分样品,并调整序列方向,条形码允许的错配数为0,最大引物错配数为2。DNA 抽提、荧光定量PCR、MiSeq 测序、序列质控、拼接均由商业测序公司完成(美吉生物,上海)。

1.6 数据分析

以97%的相似度[25]对序列进行操作分类单元(OTU)聚类并剔除嵌合体,以RDP classifier比对Silva 16S rRNA、18S rRNA数据库(v138.1)来对每条序列进行物种分类注释,设置比对阈值为70%。获得OTU数据表后,分别筛选原核生物域和真菌界,分别筛得24380、1493个OTU。保留至少3个样本中reads数≥5的OTU,并按最小样本数抽平。所得结果中,每个重复组细菌的OTU频数合计为23698,测序深度指数≥96.0%；真菌的OTU频数合计为29965,测序深度指数≥99.8%；测序深度指数表明处理得到的数据集足以代表该区域真实的微生物群落。对细菌群落和真菌群落,分别将其Sobs指数、Shannon多样性指数以及丰度占比前8的门类的各样点reads数等10个因子组合作为微生物群落组成因子(bacteria, fungi)。

细菌物种的功能属性的文献查阅以原核生物标准命名列表(LPSN)数据库的信息为主,若一个细菌属下所有公开发表(validly published)的物种均有耐盐3%及以上或生长所需盐浓度的下限>0.5%,则该属被识别为耐盐菌；若一个细菌属所有公开发表的物种均具备反硝化能力,则该属被识别为反硝化菌；反硫化菌和自养菌等鉴定方法同上。

各扇区顶级景观植被群落的差异可能体现成土时间对演替阶段的限制。本文对各扇区顶级景观植被群落的微生物进行PICRUSt2(Phylogenetic Investigation of Communities by Reconstruction of Unobserved States 2)分析,再将取得的基于KEGG(Kyoto Encyclopedia of Genes and Genomes)数据库的功能通路信息结果进行线性判别分析(LDA),识别对各样地功能通路组成的差异效果影响较大的通路,判别标准为LDA得分(log10转化)大于2,且克鲁斯卡尔-沃利斯秩和检验(Kruskal-Wallis rank sum test)为统计学显著(LDA score>2,P<0.05)。

对成土时间的驱动力进行偏曼特尔检验(Mantel)时,将本文中除成土时间外的其他环境因子(distance、vegetation、soil properties、bacteria、fungi)构建控制环境矩阵(算法为欧氏距离),检验植被、细菌、真菌群落的物种矩阵(算法为布雷柯蒂斯距离的可变性能够被成土时间解释的部分。环境变量和生物变量相关分析(Bioenv)选取本研究除bacteria、fungi之外的18个环境因子。VPA分析时,由于微生物OTU矩阵0值较多,在计算Bray-Curtis距离时进行了海林格距离转化,绘图方法如Delgado-Baquerizo等的研究[6]。

多样性统计、Bioenv分析、Mantel分析、方差分解分析使用R包vegan完成。相关性分析使用R软件包ggcor完成。LDA分析和LEfSe分析于Galaxy/Hutlab网页完成。物种Venn图和PICRUSt2在美吉生信云网页绘制。其他统计图使用ArcGIS 10.2、Adobe Illustrator 2021和R软件包ggplot2、reshape2绘制。

2 结果与分析

2.1 土壤理化性质和微生物OTU组成异质性检验

在百年尺度上,黄河三角洲流域不同时期扇区的土壤理化性质相近(图2)。河口海岸距离与有机质(r=0.75)、全氮(r=0.65)、电导率(r=-0.82)、pH(r=-0.42)极显著相关(P<0.001),且以上4个土壤理化指标相互间均有线性耦合(P<0.001)。成土时间仅与速效钾显著相关(r=-0.27,P<0.05),但非影响速效钾的首要环境因素,速效钾主要受到土壤有机质(r=0.40,P<0.01)和全氮(r=0.48,P<0.001)的影响。

图2 成土时间、海陆距离与土壤理化指标的相关性特征

各扇区微生物群落的OTU组成相近(图3)。细菌群落共筛选出6374个OTU,其中4381个(68.7%)为4个扇区共有OTU,272个(4.3%)为单个扇区特有OTU(图3)。真菌群落共筛选出391个OTU,其中260个(66.5%)为4个扇区共有OTU,15个(3.8%)为单个扇区特有OTU。1512个(23.7%)细菌OTU的丰度随成土时间的变化显著(P<0.05),67个(17.1%)真菌OTU的丰度随成土时间的变化显著(P<0.05)。 4个冲积扇区自然禀赋相似,土壤微生物OTU组成相近,生态系统的空间演替可充分体现成土时间对生物演替的驱动作用。

图3 各冲积扇区细菌和真菌OTU数量分布特征

2.2 植物分布特征分析

伴随成土时间的增加,植物物种数和科类丰富度整体呈增加趋势(表2)。研究区域内,代表性的盐生植物有芦苇(phragmitescommunisTrin.)、柽柳 (TamarixchinensisLour.) 及翅碱蓬 (SuaedaheteropterKitag)等,优势植物类群为禾本科(Gramineae)、藜科(Chenopodiaceae)、菊科(Compositae)、柽柳科(Tamaricaceae)等。中生草地群落的指示种为白茅(Imperatacylindrica(L.) Beauv.)、罗布麻(ApocynumvenetumL.)、野大豆(GlycinesojaSieb. et Zucc.)等。成土年代最近的R4扇区尚未发育出灌丛草地的植被类型,而最早的R1扇区出现了自然生长的白蜡、洋槐等树木。

R4S4、R3S5、R2S4、R1S4为各冲积扇区发现的内陆顶级景观植被(表1),建群种除柽柳外,从耐盐的狗牙根、芦苇向轻度耐盐[8]的茵陈蒿和白茅过渡(表2)。R1S3、R2S4、R3S4和R4S3为各扇区中科类丰富度较高的样地,且一致分布在距海20km附近。随成土时间增加,距海20km附近依次呈现从盐生草地、中生草地向灌丛草地或高禾草地景观的过渡特征。

表2 各样地植被物种多样性

成土时间驱动植物生物量(r=0.53,P<0.001)和植物物种丰富度(r=0.29,P<0.05)增加,且植被盖度、Shannon多样性均随成土时间增加而增大(图4)。内陆区域受人类活动干扰的影响,农田植物物种多样性较低,表现为植物Shannon多样性与土壤有机质和全氮呈负相关(P<0.01)。

2.3 细菌和真菌群落组成分布特征

研究区域内,各扇区微生物物种门类组成相似,但各门类相对优势度存在差异(表3)。细菌群落3个主要优势门类为变形菌门(Proteobacteria)、放线菌门(Actinobacteria)、绿弯菌门(Chloroflexi),总占比为59.7%。真菌群落3个主要优势门类为子囊菌门(Ascomycota)、担子菌门(Basidiomycota)、毛霉门(Mucoromycota),总占比为94.9%。

表3 研究区域内各冲积扇区细菌和真菌门类相对丰度/%

细菌群落的Shannon多样性受到土壤pH限制(r=-0.26,P<0.05,图5)。有机质含量显著促进真菌群落物种数(Sobs指数)增加(r=-0.62,P<0.001),电导率和pH则是限制真菌物种数的限制因子(r=-0.69,P<0.001,图5)。植被盖度和植物Shannon多样性的提高促进真菌物种多样性增加(r>0.28,P<0.05)。

图4 成土时间、海陆距离、土壤理化指标与植物多样性的相关性特征

图5 细菌和真菌多样性、门类相对丰度与环境因子的相关性热图

河口海岸距离、土壤有机质、全氮、pH和电导率是影响微生物门类组成的主要因子,与多数主要门类的相对丰度有较强的相关性(图5)。成土时间对微生物群落门类相对丰度影响则相对独立。

随着时间演替,部分门类的占比变化呈现相反趋势。成土时间驱动硝化螺菌门(Nitrospirota)、酸杆菌门(Acidobacteriota)、Methylomirabilota相对丰度显著增加(r>0.29,P<0.05,图5),同时拟杆菌门(Bacteroidota)、芽单胞菌门(Gemmatimonadetes)、Patescibacteria相对丰度显著降低(r<-0.27,P<0.05)。

成土时间还驱动Sumerlaeota、装甲菌门(Armatimonadota)、异常球菌门(Deinococcota)、Acetothermia、Hydrogenedentes、盐厌氧菌门(Halanaerobiaeota)、海微菌门(Marinimicrobia-SAR406-clade)、Calditrichota、Campilobacterota等次要门类相对丰度的降低(r<-0.26,P<0.05,图4未列出)。

丰度与成土时间相关性显著的真菌门类为捕虫霉门(Zoopagomycota)(r=0.30,P<0.05)。研究区域内,子囊菌门(Ascomycota)和担子菌门(Basidiomycota)在真菌群落的丰度占比达91.9%(表3)。对这2个门类的相对丰度影响最大的环境因子为pH(r=0.60,P<0.001,图5)。

2.4 细菌和真菌群落差异物种分析

在从门到属的分类层级,部分微生物群落物种类群富集特征显著(LDA score>2,P<0.05,图6)。如,细菌厚壁菌门(Firmicutes)和真菌毛霉门(Mucoromycota)和捕虫霉门(Zoopagomycota)在较早形成的R1扇区中可作为真菌群落的生物标志物(biomarker),而担子菌门(Basidiomycota)在R4扇区丰度最高(LDA score>2,P<0.05)。

图6 各冲积扇区细菌和真菌LEfSe生物标志物

富集于R3或R4扇区的细菌属,包括Erythrobacter[26]、Fusibacter[27]、Gracilimonas[28]、Halanaerobium[29]、Limibaculum[30]、Marinobacter[31]、Palleronia-Pseudomaribius[32]、Pontibacter[33]、Salinimicrobium[34]、Zeaxanthinibacter[35]等,经文献分析发现除Erythrobacter[26]外均为耐盐类群。富集于R4扇区的Campilobacterota为硫氧化化能自养细菌[36—38]。富集于R3或R4扇区的大型真菌为盘菌纲(Pezizomycetes)和烟管菌属(Bjerkandera),而富集于R1的大型真菌为隶属于蘑菇科(Agaricaceae)的侧耳属(Pleurotus)。以上物种均对各冲积扇区的微生物群落结构组成的差异具有重要影响。

为控制海陆距离变量,对距海20km附近的R1S3、R2S4、R3S4和R4S3等4个样地的细菌群落的属进行多级物种差异判别分析(Linear discriminant analysis Effect Size,LEfSe)。去除“norank”或“unclassified”开头的没有明确分类信息或分类名称的类群,并去除平均丰度为0(Mean=0)的属,最终得到554个细菌属。215个(38.8%)属对距海20km附近微生物群落结构组成的差异具有重要影响(LDA score>2,P<0.05),其中73个属的显著富集区域为R4S3样地,其他68、21、53个属分别为R3S4、R2S4和R1S3样地的生物标志物(biomarker)(表4)。在这215个属中,至少100个属为嗜盐细群,其中有86个属(86.0%)的显著富集区域为成土较晚的R3和R4扇区；涉及反硝化、反硫化、光异养、嗜热等功能性状的细菌类群同样集中于R3和R4扇区；R4扇区显著富集的自养微生物多达13种,其中10种可营化能自养。

表4 距海20km 附近215个组间丰度差异显著细菌属的部分功能属个数

2.5 各扇区顶级景观群落细菌和真菌基于PICRUSt2通路丰度的LDA判别分析

针对细菌通路, PICRUSt2分析了各扇区顶级景观群落功能通路的异质性。各样地细菌群落代谢通路相似, 但是代谢通路的丰度值不同,其中有12个level 2代谢通路存在组间显著差异(图7)。

R4S4样地显著富集的细菌群落功能通路较多涉及基础生命活动。在该样地的细菌群落,碳水化合物代谢(Carbohydrate metabolism)和氨基酸代谢(Amino acid metabolism)两个level 2代谢通路的丰富度相对较高(LDA score>2,P<0.05,图7)。存在组间显著差异的level 3代谢通路中,所有从属于遗传信息处理(level 1: Genetic Information Processing)的9个通路均富集于R4S4扇区,涉及氨酰tRNA生物合成(ko00970)、核糖体(ko03010)、RNA降解(ko03018)、DNA复制(ko03030)、蛋白质运出(ko03060)、核苷酸切除修复(ko03420)、错配修复(ko03430)、同源重组(ko03440)和内质网蛋白质加工(ko04141)等生命活动必需的复制和表达过程(LDA score>2,P<0.05,表5)。此外,涉及嘌呤代谢(ko00230)、嘧啶代谢(ko00240)、氨基酸生物合成(ko01230)等关键分子的通路的相对丰度均显著高于其他扇区(LDA score>2,P<0.05)。光合作用(ko00195)和光合生物固碳(ko00710)通路在R4S4样地细菌群落相对丰度较高(LDA score>2,P<0.05)。

R1S4样地细菌群落富集的细菌群落功能通路则涉及多样化的代谢过程。涉及细胞老化(Aging)、外源物质的生物降解与代谢(Xenobiotics biodegradation and metabolism)、复制与修复(Replication and repair)等level 2通路丰度高于其他样地(LDA score>2,p<0.05,图7)。在营养代谢方面,R1扇区有关单糖的通路丰度较高,如糖酵解/糖异生(ko00010)、戊糖和葡萄糖醛酸转化(ko00040)、果糖和甘露糖代谢(ko00051)、半乳糖代谢(ko00052)、丁酸(ko00650)等(LDA score>2,P<0.05,表5)。芳香族化合物降解(ko01220)的代谢功能也富集于R1S4样地(LDA score>2,P<0.05)。

表5 细菌群落从属基因信息处理和4类代谢通路的KEGG level 3差异通路

针对真菌通路, PICRUSt2分析得到了各扇区顶级景观群落MetaCyc pathway信息及其丰度,共有25条MetaCyc代谢通路存在组间显著差异(图7)。R4S4样地显著富集的真菌代谢通路亦较多涉及基础代谢,如TCA循环(PWY-5690)、丙二醇生物合成(PWY-7385)、有关细胞色素c的有氧呼吸I(PWY-3781、PWY-7279)和葡萄糖酵解(ANAGLYCOLYSIS-PWY)通路(LDA score>2,p<0.05,图7)。R1S4样地显著富集的主要是有氧呼吸相关的泛醌类合成通路,如泛醌- 6(PWY-7235合成)、泛醌- 7(PWY-7235)和泛醌- 9合成(PWY-7235)(LDA score>2,P<0.05)。

图7 细菌和真菌群落KEGG途径丰度差异的线性判别分析(LDA)条形图

2.6 成土时间驱动力的统计学检验

成土时间对植物、细菌、真菌的群落组成变异的解释能力分别为4.76%、2.98%、3.68%(图8),均大于全球范围土壤年龄对生态系统指标的平均解释能力2.1%[6],表明成土时间在黄河三角洲是重要的生态驱动力。成土时间对植物、细菌、真菌的群落组成变异的独特方差解释比例分别为3.41%、2.61%、1.12%,独特解释比例大于0,说明成土时间对生物群落组成变化的影响不能被其他状态因子完全代替,因此成土时间是为黄河三角洲生物演替带来额外解释度的显著驱动因子。

图8 植被和微生物群落组成的方差分解分析

成土时间对细菌和真菌群落组成矩阵均显著偏相关(Partial MantelP<0.05)(表6),表明成土时间单一因子对微生物演替过程中群落物种组成有极显著驱动力。若不加入成土时间,细菌和真菌群落的Bioenv模型的相关系数相对最佳模型的拟合精度分别相对降低1.41%、5.78%(表7),表明在黄河三角洲微生物演替研究中将成土时间纳入考虑可提高拟合精度。

3 讨论

3.1 成土时间对黄河三角洲植物群落演替的影响

由于成土时间带来演替时间的梯度,各扇区内陆植物群落结构呈由简单向复杂的演变趋势。研究区域内,植物分布和多样性在随成土时间的差异主要体现在生物量(r=0.53,P<0.001)和植物物种数(r=0.29,P<0.05)。距海约20km的自然植被样地呈现了植被类型由低级向高级的过渡特征。生态系统的空间演替体现了成土时间对自然植被由低级向高级发展的驱动作用。刘志杰[39]的样方调查同样发现,成土时间在1964年后的黄河三角洲叶瓣植物种类较少,植物种类不超过10科；灌草丛在成陆时期较长的区域成为主要植被类型之一。陈彭祯霓等[40]以冰川退缩区百年时间跨度的植被初生演替为研究对象,发现不同演替时间的样地植物物种组成差异较大,木本植物生物量随演替时间增加而增加；黄河三角洲植物群落结构演替规律与之相符,丰富了百年尺度的植被初生演替规律的研究。

随成土历史的延长,白茅的分布向海岸推进,体现了成土时间对近海土壤生境条件改善的驱动作用。白茅的出现代表演替到达较高阶段,种群结构趋于稳定[41]。本研究中,在白茅分布地点,植被的共同点为植被盖度在70%以上。在R1扇区,白茅为常见种。在R4、R3和R2冲积扇区,白茅分别在距入海口约20km、14km及9km进入群落并成为伴生种。R4冲积扇区未演替出白茅占优的草地群落,白茅仅在人类聚落附近偶见；R3S4样地是在R3扇区发现的距海最近的白茅群落,白茅盖度约25%；而在距海仅15.8km的R2S3样地,白茅盖度已达40%。前人研究证实,白茅是黄河三角洲盐土改良的指示植物,如侯本栋等[2]发现,白茅群落的土壤有机质、全氮含量显著高于柽柳群落、盐地碱蓬群落和裸地,且含盐量相对较低。

3.2 成土时间对黄河三角洲微生物演替的影响

成土时间的增加可驱动细菌和真菌群落结构(OTU的Bray-Curtis距离矩阵)的改变,这种改变具有一定的方向性(Partial Mantelr≥0.1377,P<0.05,表6)。除去土壤理化性质、植物互作关系等的共线性干扰后,成土时间对微生物群落变化依然具有独特影响；成土时间可为群落变化趋势提供额外的解释度(图8,表7)。成土时间在黄河三角洲局地尺度是独特而显著的重要生态驱动力(图8),这与Delgado-Baquerizo等[6]对土壤年龄驱动力的整合分析相符,即成土时间主要影响与生物活动相关的缓慢变化的生态指标。

表6 成土时间对生物群落组成变化影响的偏Mantel检验

表7 成土时间对微生物群落Bioenv最佳环境因子组合的拟合系数的影响检验

微生物群落的功能特征也随成土时间增加而发生显著改变(图7,表5)。在成土时期较近的区域,微生物群落涉及基因信息处理等基础生命活动的功能通路相对丰度较高,与Yang等[42]关于黄河三角洲氨基酸代谢通路在高盐土壤丰度更高的发现一致,推测在演替初期环境胁迫下,微生物群落依靠功能冗余提高抗逆性[43]及维持群落功能的稳定性[44—45]；也可能是响应于频繁而幅度较大的环境变化,微生物需要能够迅速利用周边资源生存、生长和繁殖[21],故采取r-选择竞争策略[46]。在生态系统年龄较长的R1扇区,涉及细胞老化调控等功能通路的丰度显著较高(LDA score>2,P<0.05)符合环境趋于稳定带来选择过程优势的增加[47]的理论；单糖代谢、芳香族化合物降解等碳代谢相关基因通路的丰度发生显著变化(表5),为前人关于微环境因素影响微生物碳源利用[42, 48]、土壤发育对生态位分化有促进作用[21]等论点提供了支撑。

功能微生物的富集特征是微生物群落组成随成土时间变化的原因之一。LEfSe分析表明,在各扇区距海20km附近的各样地,对细菌群落组成变异贡献显著的215个可识别属中,嗜盐细菌、反硝化细菌、反硫化细菌等功能微生物的显著富集区域均为成土较晚的R3和R4扇区(表4),其优势度随演替时间有下降趋势。对不同扇区细菌群落组成的整体差异起关键作用的细菌属也以耐盐菌为主,集中区域同样为成土较晚的R3和R4扇区(图6)。范运梁[49]依据微生物培养实验发现,伴随盐生植物的生长,耐盐细菌的比例显著下降。苗菁等[11]还发现,随着植被的覆盖以及植物源有机物的输入,盐度对土壤细菌属的作用逐渐减弱。耐盐细菌类群的减少趋势与盐度随演替降低的趋势[8, 21]相适应。随盐土向农田土壤演替,根系生长可改善土壤结构、降低土壤容重、增加土壤孔隙[50],通气性的增加促进好氧型生物活性的增加[51],而反硝化、反硫化均为厌氧过程[52—53],导致反硝化菌与反硫化菌竞争力随成土时间增加的下降。

成土时间引起部分适应高营养环境的类群相对丰度的增加,以及部分先锋类群相对丰度的降低。酸杆菌门(Acidobacteriota)在环境土壤中一般占细菌总量20%左右[54],在土壤营养较高的黑土地或森林土壤,酸杆菌门占比更高[55—57]。酸杆菌相对丰度平均占比仅为8.26%(表3),苗永君[12]在黄河三角洲细菌多样性的研究同样发现,酸杆菌在农用地和自然植被土地类型的相对丰度占比均<10%。而且,酸杆菌随土壤有机质增加(r=-0.47,P<0.001)和电导率下降(r=-0.39,p<0.01)而增加(图5),表明其在黄河三角洲属于随生境改良而增加的类群。拟杆菌门(Bacteroidota)、Sumerlaeota、Hydrogenedentes、海微菌门SAR406分支(MarinimicrobiaSAR406-clade)、异常球菌门(Deinococcota)和浮霉菌门(Gemmatimonadota)则随土壤有机质增加(r<-0.33,P<0.01)和电导率降低(r>0.33,P<0.01)而出现相对丰度的显著下降,是优势度随土壤改良而降低的先锋类群。在成土早期相对贫瘠的极端生境,适应寡营养生境的先锋微生物容易取得较高优势度[10, 21],但是土壤营养限制的解除不利于这些种群的壮大[58],如李辉等[59]发现,在赤泥自然成土的过程中,富营养细菌与贫营养细菌丰度比值随演替明显增加。

成土时间的增加驱动微生物生态位的分化。细菌群落功能通路的分化趋势(图6)体现了微生物生态位的分化。营养方式的多样化同样是微生物生态位分化的体现,如寄生土壤线虫、轮虫和变形虫的捕虫霉门[60](Zoopagomycota)真菌(r=0.30,P<0.05),寄生昆虫的黑僵菌[61](Metarhiziumanisopliae)真菌(r=0.32,P<0.05),寄生变形虫的Neochlamydia[62]细菌(r=0.26,P<0.05),与共生固氮的Ensifer[63](r=0.31,P<0.05)和Azoarcus[64](r=0.33,P<0.01)细菌均为优势度随成土时间增加而增大的类群,演替初期此类物种受制于宿主资源的贫乏。在对距海20km附近样地细菌群落组成的对比中,还发现与反刍动物共生的瘤胃球菌属[65](Ruminococcus)在R1S3样地的显著富集(LDA score=2.28,P<0.05),说明随演替时间而增强[3]的人类活动能为微生物提供更多可选生态位。

3.3 成土时间驱动力的微生物因素

植被是滨海生态系统发展其生态系统功能最关键的过程[10]。盐生植物盖度提高可减弱土壤的自然蒸腾,抑制盐分在土壤表层的持续累积；翅碱蓬等真盐生植物能大量吸收盐分,使土壤脱盐[66]。盐生植被改善土壤质量的作用为微生物提供良好的生存条件,植物根系的分泌物及植物残体是土壤微生物的主要可利用的碳源[58, 67]。

但是,成土初期盐滩光板地的存在表明先锋植物的定居需要等待土壤的形成和改良过程。在极端生境下,高等植物无法存活,但细菌群落物种数和多样性较少受环境因子影响(图5,表明某些适应力较强的菌种可定居。在微生物定居于原生裸地的同时,即开始了成土进程[68],成土初期的功能微生物促进了生境的改良。

自养微生物促进土壤碳氮积累。光合作用相关通路富集在成土相对较早的扇区,这与Wang等[10]光合作用相关基因显著富集于新生湿地的发现一致。蓝藻细菌门能够进行产氧的光合作用, 部分物种能进行生物固氮,还可以通过光合作用释放特异的蛋白质信号分子来调控植物的生长和发育[67]。但是,本研究中蓝藻门(Cyanobacteria)相对丰度与成土时间无显著相关性(P>0.05,图4),暗示蓝藻以外的光合生物贡献于细菌群落光合作用表型的富集。本研究在16S rDNA克隆文库中发现了CandidatusChloroploca[69]、Rhodobacter[70]等光合硫氧化细菌以及光异养的Rubribacterium[71]、Rubellimicrobium[71]等紫色非硫细菌,还发现了距海20km附近发现的光异养细菌在R3、R4扇区的富集(表4),此类光合生物很可能对成土早期的细菌光合作用起到了关键作用。在R4S3扇区中还发现了化能自养细菌的富集现象(表4),包括3种耐盐的自养硝化菌和6个可硫氧化自养的细菌属。前人强调成土早期以蓝藻为主的光合自养微生物对土壤碳氮的积累起到关键作用[59, 72],本研究推测滨海盐土化能自养菌和蓝藻以外的光合生物在土壤碳氮积累方面的作用同样具有生态学意义,这与近海河流沉积物硫元素的富集和活跃的元素化学循环[73]相适应,有待实验验证。除蓝藻和自养硝化细菌外,耐盐固氮菌也具有促进土壤氮积累的功能潜力。16S rDNA克隆文库中发现的耐盐的共生固氮菌Ensifer[63]和Azoarcus[64]分别可以与豆科植物和禾本科植物形成根瘤；耐碱的中慢生根瘤菌属[74—75](Mesorhizobium)也有为先锋植物提供了早期土壤所缺乏的氮源的潜力。

李辉等[59]还发现,微生物可通过分泌有机酸、菌丝改良土壤结构等方式参与土壤的改良。微生物能够释放土壤酶, 催化氧化还原、有机质矿化、腐殖物质在土壤中的合成、以及生长活性物的释放等多种生物化学反应[22], 改善植物的生长环境。盐生植被与微生物循环反复作用[18], 最终实现滨海盐土向农田的转化。

4 结论

成土时间对细菌和真菌群落变化具有显著的偏相关性(MantelP<0.05),对植物、细菌和真菌群落组成变化的解释能力高于全球尺度土壤年龄对生态系统指标的平均解释能力。成土时间驱动黄河三角洲植物、细菌和真菌群落组成变化,也造成黄河三角洲微生物群落的功能特征的显著异质性。在百年尺度上,成土时间是黄河三角洲植物、细菌和真菌群落演替的重要而显著的驱动力,将成土时间纳入研究变量能够为生物群落演替的研究提供额外的解释度。该结论为提高评估滨海生态系统植物和微生物演替规律的准确性提供了科学依据。