Types and diversity of natural regeneration community after clear cutting of Phyllostachys edulis forests in Mount Tianmu, China

浙江天目山国家级自然保护区位于浙江省西北部杭州市临安区境内的西天目山(30°18′30″~30°24′55″N，119°23′47″~119°28′27″E)，总面积为1 050 hm²，主峰海拔为1 506 m。年平均气温为8.8~14.8 ℃，年积温为2 500~5 100 ℃，年降水量为1 390~1 870 mm，相对湿度为76%~81%^[7]。自然保护区受海洋暖湿气候影响较深，森林植被十分茂盛。由于区内地势较为陡峭，海拔上升快，气候差异大，植被的分布有着明显的垂直界限，自山麓到山顶垂直带谱为：海拔870 m以下为常绿阔叶林；870~1 100 m为常绿、落叶阔叶混交林；1 100~1 380 m为落叶阔叶林；1 380~1 506 m为落叶矮林^[8]。自然保护区中的毛竹林主要分布于海拔350~900 m其他森林类型之中，毛竹林下植被稀少，有少量的豹皮樟Litsea coreana var. sinensis、短尾柯Lithocarpus brevicaudatus、微毛柃Eurya hebeclados、细叶青冈Cyclobalanopsis myrsinaefolia、连蕊茶Camellia fraterna、马银花Rhododendron ovatum和牛鼻栓Fortunearia sinensis等常绿树种^[9]。

2018年11−12月，在浙江天目山国家级自然保护区内2015年底的毛竹林皆伐试验区内全面踏查的基础上，采用样方法，设置11个20 m × 20 m群落样地，在每个群落样方四角及中间位置分别设置1个5 m × 5 m木本层样方及草本样方，共55个木本层样方和55个草本样方。

对样方内的木本植物进行每木调查，记录每个植株的物种名、胸径、树高、冠幅等并挂牌以便长期监测。记录每个草本样方内草本层植物的物种名、多度、盖度、平均高度、平均基径。并用全球定位系统(GPS)和罗盘仪测量每个样方的经度、纬度、海拔、坡度、坡向等环境因子(表1)。

物种重要值是表示物种在群落中优势度的一个指标，可以表示不同物种在群落中的功能地位^[10]。分别计算每块样地的重要值，计算公式为^[11-12]：I_V灌=(相对高度+相对盖度+相对密度)/3；I_V草=(相对高度+相对盖度+相对密度)/3；相对密度=(某个种的个体数/全部种的个体数)×100；相对高度=(某个种的平均高度/所有种平均高度之和)×100；相对盖度=(某个种的盖度/所有种盖度之和)×100。

利用R语言Vegan包中的diversity函数计算各层次α多样性指数，公式如下：

Shannon-Wiener 指数H^[13]：H = −$\sum\limits_{i = 1}^S {} $P_i ln P_i；Simpson 多样性指数D^[14]：D = 1−$\sum\limits_{i = 1}^S {} $P_i²；Gleason 丰富度指数R：R=S/lnA；基于Shannon-Wiener 指数的Pielou 均匀度指数：E=H/lnS。其中，S为样方内物种数，P_i为种i的相对重要值，A为样方面积。采用R3.5.3软件中的Vegan包计算各个样方α多样性，采用Excel 2003软件对多样性数据进行统计处理，采用SPSS 25统计分析软件中的one-way ANOVA(单因素方差分析) 分析不同群落之间的α多样性的差异，差异程度用Ducan的多重比较SSR(P＜0.05)。

TWINSPAN分类运用计算机软件TWINSPAN for Windows V2.3完成^[15]。为了消除稀有种对分类结果的影响，在对天目山毛竹林皆伐后次生群落进行TWINSPAN分类时将物种在样方内出现次数≤4且所在样方重要值小于等于0.1的物种排除。最终用于TWINSPAN分类的物种数由318种删减至116种，最后将这116种植物的重要值数据与55个样方共同组成植被矩阵进行分类。

多样性指数排序使用R3.5.3软件中的Vegan包完成，排序图绘制使用ggplot2包完成。在排序前需要将测得的0~360°的罗盘仪测量值使用TRASP法转换为0~1的值。转换公式如下^[16]：T_RASP={1−cos[(π/180)(a_spect−30)]}/2。其中：T_RASP为坡向转换值；a_spect为坡向方位角角度。经过转换后坡向转换值的范围为0~1，转换后的值越大代表的样方生境越向阳。选择排序方法前首先对多样性数据矩阵进行去趋势对应分析(DCA)，如果DCA排序前4个轴中最大值超过4，选择单峰模型更合适。如果介于3~4，单峰模型和线性模型都可行。如果是小于3，则选择线性模型。经计算在DCA分析中木本层和草本层的最长排序轴长度分别为0.435 49和0.660 27，结果均小于3因此选择线性模型中的冗余分析(RDA)比较合适。因此，本研究采用RDA排序的方法将选取的地形因子作为解释变量以各多样性指数作为响应变量进行分析。

将毛竹林皆伐自然更新形成的55个样方进行TWINSPAN分类，命名原则参照《中国植被》^[17]中的群落命名原则，结合样方的实际情况，将划分的样方适当进行调整，最终划分了11个群落类型(图1)。

Figure 1.  Two-way indicator tree classification of natural restoration community after clear cutting of Moso bamboo forest in Mount Tianmu

臭椿Ailanthus altissima+掌叶复盆子Rubus chingii-栗褐薹草Carex brunnea+金星蕨Parathelypteris glanduligera群丛(Ⅰ)，包括样方46、47、48、49、50，分布于海拔829~849 m处。木本层优势种臭椿平均高度为1.8 m，盖度15%~35%；掌叶复盆子平均高度为1.8 m，盖度2%~16%；主要伴生种为檫木Sassafras tzumu、朵花椒Zanthoxylum molle、蓬蘽Rubus hirsutus等。草本层优势种栗褐薹草平均高度为0.3 m，盖度10%~20%；金星蕨平均高度为0.3 m，盖度1%~5%；主要伴生种为五节芒Miscanthus floridulus、博落回Macleaya cordata、显子草Phaenosperma globosa等。

棘茎楤木Aralia echinocaulis+朵花椒-栗褐薹草+显子草群丛(Ⅱ)，包括样方51、52、53、54、55，分布于海拔849~855 m处。木本层优势种棘茎楤木平均高度为2.2 m，盖度6%~28%；朵花椒平均高度为2.8 m，盖度15%~28%；主要伴生种为山鸡椒Litsea cubeba、檫木、掌叶复盆子等。草本层优势种栗褐薹草平均高度为0.3 m，盖度3%~8%；显子草平均高度为0.7 m，盖度5%~8%；其中栗褐薹草处于主要地位，主要伴生种为博落回、求米草Oplismenus undulatifolius、密叶薹草Carex maubertiana等。

朵花椒-博落回+显子草群丛(Ⅲ)，包括样方21、25、31、32、33、34、35，分布于海拔527~570 m处。木本层优势种朵花椒平均高度为3.4 m，盖度10%~43%；主要伴生种有臭辣树Evodia fargesii、白背叶Mallotus apelta、醉鱼草Buddleja lindleyana等。草本层优势种博落回平均高度为2.0 m，盖度3%~10%；显子草平均高度为0.6 m，盖度2%~10%；主要伴生种有栗褐薹草、庐山香科科Teucrium pernyi等。

朵花椒+白背叶-栗褐薹草群丛(Ⅳ)，包括样方22、26、27、28、29、30、37，分布于海拔430~435 m、559~596 m处。木本层优势种朵花椒平均高度为4.2 m，盖度10%~22%；白背叶平均高度为2.5 m，盖度4%~43%；主要伴生种为臭辣树以及少量掌叶复盆子、苎麻Boehmeria nivea等。草本层优势种栗褐薹草平均高度为0.4 m，盖度3%~21%。主要伴生种为显子草和博落回等。

白背叶-显子草+五节芒群丛(Ⅴ)，包括样方1、2、3、4、5，分布于海拔535~545 m处。木本层优势种白背叶平均高度2.4 m，盖度15%~27%。主要伴生种为老鸦糊Callicarpa giraldii、朵花椒、掌叶复盆子等。草本层优势种显子草平均高度0.8 m，盖度5%~15%；五节芒平均高度1.2 m，盖度3%~10%；主要伴生种有博落回、栗褐薹草、山莓Rubus corchorifolius等。

白背叶-苦竹Pleioblastus amarus+三穗薹草Carex tristachya群丛(Ⅵ)，包括样方11、12、13、14、15，分布于海拔568~578 m处。木本层优势种白背叶平均高度2.5 m，盖度15%~45%。主要伴生种为老鸦糊、浙江柿Diospyros glaucifolia等。草本层优势种为苦竹平均高度0.9 m，盖度5%~15%；三穗薹草平均高度0.4 m，盖度2%~6%；主要伴生种为栗褐薹草等。

朵花椒-青绿薹草Carex breviculmis群丛(Ⅶ)，包括样方23、24、42、44，分布于海拔559~569 m、866~876 m处。木本层优势种朵花椒平均高为4.0 m，盖度13%~58%。主要伴生种为蓬蘽等。草本层青绿薹草为绝对优势种，平均高度0.4 m，盖度5%~35%。主要伴生种为显子草、密叶薹草等。

蓬蘽+臭椿-显子草+蕨Pteridium aquilinum群丛(Ⅷ)，包括样方6、9、41、43、45，分布于海拔538~547 m、866~874 m处。木本层优势种蓬蘽平均高0.34 m，盖度4%~37%；臭椿平均高1.9 m，盖度1%~26%；主要伴生种为臭辣树、苎麻等。草本层优势种显子草平均高度为0.6 m，盖度3%~8%；蕨平均高度为0.7 m，盖度2%~12%；主要伴生种为密叶薹草、三穗薹草等。

糙叶树Aphananthe aspera+白背叶-三穗薹草+博落回群丛(Ⅸ)，包括样方36、38、39、40，分布于海拔430~435 m处。木本层优势种糙叶树平均树高4.3 m，盖度32%~52%；白背叶平均树高3.2 m，盖度5%~35%；主要伴生种为茶Camellia sinensis、高粱泡Rubus lambertianus等。草本层优势种三穗薹草平均高度为0.2 m，盖度3%~35%；博落回平均高度为2.0 m，盖度1%~10%；主要伴生种为青绿薹草、显子草等。

蓬蘽-皱叶狗尾草Setaria plicata群丛(X)，包括样方16、17、18、19、20，分布于海拔545~550 m处。木本层优势种蓬蘽平均高为0.5 m，盖度2%~52%。主要伴生种为山莓、小构树Broussonetia kazinoki。草本层优势种皱叶狗尾草平均高度0.6 m，盖度23%~37%，并伴生少量求米草、黑足鳞毛蕨Dryopteris fuscipes等。

紫麻Oreocnide frutescens+牛鼻栓Fortunearia sinensis-野茼蒿Crassocephalum crepidioides+显子草群丛(Ⅸ)，包括样方7、8、10，分布于海拔538~547 m处。木本层优势种紫麻平均高度为1.5 m，盖度58%；牛鼻栓平均高度3.6 m，盖度40%；主要伴生种为小齿钻地风Schizophragma integrifolium、木莓Rubus swinhoei等。草本层优势种野茼蒿平均高度0.5 m，盖度2%~12%；显子草平均高度0.7 m，盖度4%~8%；主要伴生种为金灯藤Cuscuta japonica等。

通过TWINSPAN分析将55个样方分成了11类群丛，其中群丛Ⅰ、Ⅱ、Ⅴ、Ⅵ、Ⅹ分类结果与样地分布一致，其他群丛在分类时出现了同属一个样地而分类结果不同的情况。TWINSPAN分析是基于群落各物种重要值的大小而进行的，因此同属一个样地但样方内植被分布不同是造成这种情况的根本原因。地形上的异质生境会对植被的分布格局产生影响^[18]，地形通过形态(起伏等)变化、经过不同地貌过程(堆积、搬运等)影响了光、热、水和土壤理化性质等因素在空间上的再分配，形成异质生境^[19]。地形单位又可以进一步划分为上部坡面和下部坡面两部分，上下部坡面又可以划分为顶坡、上部边坡、谷头凹地、麓坡、泛滥性阶地、等微地形单位^[20]。各微地形单元之间受不同地貌过程、水文、土壤形成过程、水分等环境因子存在明显差异^[21]，加上不同类型、程度干扰在微地形上形成具有不同优势种的植被分异格局^[22]。此外，除了微地形单元对生境的影响外样地中还有其他未被砍伐的高大乔木随机分布，经统计，样地内D_BH≥10 cm的乔木数量为27株，其中样方内13株。这些高大乔木会影响样方内的光照、水分等条件，导致同一样地内小尺度上生境的不同，在一定程度上影响植被的分异。

由表2可以看出：在毛竹林皆伐后自然恢复形成的不同群落之间α多样性具有显著差异。群落Ⅶ、Ⅹ、Ⅺ之间木本层的Shannon-Wiener多样性指数差异不显著，但显著低于其他群落；群落Ⅰ、Ⅱ、Ⅶ、Ⅹ、Ⅺ之间草本层Shannon-Wiener多样性差异不显著，但显著低于其他群落。木本层Simpson多样性指数与Shannon-Wiener多样性指数结果类似，表现为群落Ⅶ、Ⅹ、XI显著低于其他群落，但群落VII与其他群落差异不显著；草本层Simpson多样性指数表现为群落Ⅱ、Ⅶ、Ⅹ显著低于其他群落。木本层Pielou均匀度最大值为群落Ⅰ，其值为0.888，最低为群落Ⅺ，其值为0.676；草本层Pielou均匀度最大值为群落Ⅲ，其值为0.856。群落Ⅶ的值最低。木本层和草本层的Pielou均匀度指数并无显著关联，如群落Ⅺ木本层Pielou均匀度最低但草本层Pielou均匀度最高，群落Ⅶ木本层和草本层均为最低，可能是由于群落Ⅶ的朵花椒和青绿薹草作为优势种长势较好，最终影响了其他物种的生长而导致的。木本层Gleason丰富度指数最大为群落III其值为11.406，其次是群落Ⅵ和群落Ⅴ，群落Ⅺ的Pielou均匀度指数其值为4.453。草本层Gleason丰富度最大为群落Ⅵ其值为5.033，其次是群落VIII和群落Ⅹ，群落Ⅺ的Gleason丰富度指数最低，其值为2.175。

各群落之间不同层次之间的差异性可以看出各个群落之间的α多样性有较大差别，这可能是由于不同群落之间不同的环境因子所导致的，如群落Ⅰ、Ⅱ的海拔较高但是各层次的α多样性较低，显示了海拔对群落α多样性的影响。

对木本层和草本层α多样性进行RDA排序结果显示木本层中3个排序轴对多样性变化的总解释率为38.61%，其中排序轴1的解释率为38.2%，草本层中3个排序轴对多样性变化的总解释率为17.43%，其中排序轴1的解释率为17.02%，在木本层和草本层中排序轴1的解释率比较大，因此主要对排序轴1进行分析。采用permutest函数对RDA分析结果进行蒙特卡罗置换检验(Monte Carlo permutation test)置换循环999次后显示地形因子对木本层和草本层的α多样性分布的解释量具有显著性(P＜0.05)说明了排序效果理想。用envfit函数检验每个地形因子的显著性以及地形因子与排序轴的相关性结果见表3。

在表3中地形因子与2条排序轴对应的值代表地形因子与不同排序轴的相关性，可以看出木本层的排序轴1主要与海拔呈极显著正相关(P＜0.01)，与坡向和坡度呈不显著负相关。说明沿着排序轴1自左向右海拔越高坡度越小坡向由阳坡转为阴坡。草本层的排序轴1与海拔呈显著正相关(P＜0.05)、与坡向、坡度均呈正相关，说明沿着排序轴1向右海拔越高、坡向由阴坡转为阳坡、坡度升高。r²表示该地形因子对α多样性分布的决定系数，r²越小则代表该地形因子对多样性分布的影响越小。P_r代表显著性检验，可以看出海拔对木本层和草本层的α多样性分布有着显著影响(P＜0.05)。

由图2可得：在3个地形因子中代表海拔的线段最长，代表坡度的线段最短。说明在地形因子中海拔对木本层α多样性分布影响较大而坡度的影响最小。海拔箭头方向与排序轴1方向呈锐角，说明海拔与排序轴呈正相关，这与表3的结果相同。α多样性中Shannon-Wiener多样性指数、Simpson多样性指数、Pielou均匀度箭头在第2象限，且与排序轴1呈钝角，说明这3个指标与海拔呈负相关而与坡度呈正相关，即随着海拔的升高和坡度的降低木本层多样性降低，物种分布也越加不均匀。Gleason丰富度箭头所在位置为第3象限箭头方向，与排序轴方向呈钝角为负相关，说明木本层Gleason丰富度指数与海拔呈负相关，而与坡度呈正相关，木本层Gleason丰富度高的区域往往海拔低且坡度大。

Figure 2.  RDA sequencing of woody layer α diversity

由图3可以看出：海拔的箭头处于第1象限且线段长度最长，说明海拔是影响草本层α多样性分布的主要因子，海拔的箭头方向与排序轴1方向为锐角，说明呈海拔与排序轴1呈正相关；坡向的箭头处于第4象限且与排序轴1呈锐角，说明坡向也与排序轴1呈正相关；坡度的箭头位于第1象限但是线段很短几乎看不见，说明坡度对草本层α多样性分布影响极小。草本层总体分布情况与木本层类似，Shannon-Wiener多样性指数、Simpson多样性指数、Pielou均匀度的箭头处于第2象限且与排序轴呈钝角，说明随着海拔的升高和坡向由阴转阳，草本层的Shannon-Wiener多样性指数、Simpson多样性指数和Pielou均匀度均下降，即草本层的多样性下降，物种分布均匀程度下降。Gleason丰富度的箭头处于第3象限且与排序轴1方向呈钝角，说明随着海拔的增加和坡向转阳草本层的Gleason丰富度也呈下降趋势。

Figure 3.  RDA sequencing of herbaceous layer α diversity

运用TWINSPAN方法分析结合野外调查结果及保护区中不同植物群落特征，将天目山毛竹林皆伐后自然更新形成的次生群落通过55个样方数据划分为11个群落类型，推测由于样地内微地形不同和原有高大乔木分布的不均匀形成的异质生境导致植被分类的分异。毛竹林皆伐后自然恢复形成的不同群落之间α多样性具有显著差异，具体表现为：木本层Shannon-Wiener多样性群落Ⅶ、Ⅹ、Ⅺ之间差异不显著但显著低于其他群落，草本层Shannon-Wiener多样性群落Ⅰ、Ⅱ、Ⅶ、Ⅹ、Ⅺ之间差异不显著但显著低于其他群落，这可能是由于样地之间海拔的梯度差异导致了α多样性的差异。通过分析也可以看出：海拔是影响α多样性的主要因子并与其呈负相关。

群落构建规则认为：环境筛选和生物相互作用在不同尺度上决定了群落结构、组成和多样性^[23-24]。植物群落的有序更替过程被称为演替，它以优势种的更替作为表征，环境的演变为基础，通常是一个植物群落代替另一个植物群落的过程^[25]。在演替过程中，与前期演替出现的物种相比，后期演替出现的物种具有低增长率、长寿命和耐荫等特征^[26]。本研究的群落属于演替初期，在样地内木本层共调查到落叶树种4 639株，常绿树种2 684株，木本层常绿树种占比为36.65%，且大多为蓬蘽、茶、寒莓Rubus buergeri等灌木，少部分能够在未来形成乔木层的常绿树种有小叶青冈Cyclobalanopsis myrsinifolia、杉木Cunninghamia lanceolata等，但是这部分树种占比不高，所以预计在几年后乔木层出现后会形成一个落叶阔叶林，优势种为白背叶、朵花椒一类的阳性物种。演替早期的物种大多具备种子质量小，依靠风力或鸟类传播，具有长距离扩散的能力，也可以在土壤中长期休眠，遇到地表扰动则大量萌发^[27-28]。目前，很少有对于花椒种子传播方式以及是否具有休眠特性的研究，国内对与朵花椒同属植物椿叶花椒Zanthoxylum ailanthoides的种子研究显示：椿叶花椒种子具有生理休眠的特征^[29]，而白背叶在荒地经常作为先锋树种。推测由于皆伐造成的地表扰动导致了生境的变化，使得白背叶和朵花椒可以快速占据优势位得到快速发展。在演替过程中常绿树种分布于木本层下层，在足够长的时间之后形成天目山自然保护区在这个海拔高度的地带性常绿阔叶林。

环境因子与生物多样性格局的关系一直以来都是生态学家感兴趣的领域^[30]。本研究把海拔、坡度、坡向3个地形因子作为解释变量，将自然更新群落各层次的α多样性作为响应变量，通过冗余分析(RDA)探究这两者之间的关系。经蒙特卡罗置换检验后，解释变量对木本层和草本层的α多样性相应变量都具有显著影响(P＜0.05)，因此结果可信。海拔是特殊的环境因子，反映了温度、光照和水分的综合情况^[31]。本研究通过排序轴与环境因子的关系可以看出：海拔与木本层α多样性呈极显著负相关(P＜0.01)，与草本层呈显著负相关(P＜0.05)，坡向和坡度与α多样性的各类测度指数的相关性均不显著。本研究显示更新的次生群落的木本层和草本层的α多样性都随着海拔的升高而降低这与大多数研究的结果一致^[32-33]。也有研究显示：随着海拔的升高α多样性先升高后减小^[34]或者随着海拔升高而增大^[35]。木本层排序轴解释率高于草本层的解释率的原因可能是由于草本层α多样性与土壤有机质、土壤碱解氮大气温度等微域或局部环境性状关系密切^[36]，而对地形因子变化响应不明显导致的。群落的多样性是由许多环境因子共同作用的，如pH、有机质和全氮含量^[37]等。目前，群落环境因子对多样性影响的解释率不够高，说明仍有许多环境因子没有被发掘，而在不同森林恢复阶段驱动物种多样性的环境因子也是不同的^[38]，因此，未来将进一步探索更多的环境因子对群落不同恢复阶段各层次多样性的影响。

综上所述，本研究以天目山国家级自然保护区毛竹林皆伐后自然更新群落为研究对象，采用TWINSPAN分类结合RDA排序的研究方法，取得了良好的分类结果，并揭示了演替初期群落与环境因子的分布格局。可为天目山国家级自然保护区植被恢复过程研究提供理论依据。