广东南岭20 hm2样地华南五针松空间分布与生境特征研究

龚粤宁， 刘志发， 水坤春， 张 强， 吴林芳， 郭腾辉， 权 擎， 沈 勇*

( 1. 广东南岭国家级自然保护区管理局， 广东 韶关 512727； 2. 广东省科学院动物研究所， 广东省动物保护与资源利用重点实验室， 广东省野生动物保护与利用公共实验室， 广州 510260； 3. 广州林芳生态科技有限公司， 广州 510520 )

华南五针松()，又名广东松，属松科(Pinaceae)松属()常绿乔木。2021年国家重点保护野生植物名录将华南五针松列为中国二级保护植物，其种群在湖南、广东、广西、贵州、海南等地均有分布(沈燕等，2016)。迄今为止，国内发现的面积较大、保存较完好的华南五针松种群位于广东南岭国家级自然保护区内石坑尾(112°56′8″—113°4′18″ E、24°30′28″—24°48′9″ N)，该区域分布约1 300 hm华南五针松原生林(张璐等，2006)。华南五针松目前处于近危状态(张璐等，2006)，加强该物种的研究对于物种的保护具有现实意义。目前，针对该物种的研究涉及种间关系(丁忠江等，2000)、群落特征(王献溥和李信贤，1989；古炎坤等，1993；沈燕等，2016)、种群现状(杜道林等，1996)、生理特征(段小平和陈卫军，1992；周佑勋和段小平，1993)及对气候变化的响应(黄蕴凯等，2021)等方面。因此，针对珍稀濒危植物的保护探索其空间分布格局，从而确定其群落结构形成的生态学过程(Boyden et al., 2005)，并明确其生境特征，对于制定相应的就地保护方案，以及完善其迁地保护策略具有重要的指导作用。然而，目前尚未见有针对华南五针松开展此类研究的报道。

空间点格局分析(spatial point pattern analysis)被认为是研究群落结构形成机制的有效手段，对揭示植物群落的动态具有重要作用(de la Cruz et al., 2008)。比如，空间聚集假说(the spatial segregation hypothesis)认为种内的聚集导致种间的隔离，从而导致竞争主要发生在种内的个体之间，而不是发生在种间不同个体之间。该假说预测，经过该过程，竞争能力过强的物种将会被抑制，避免处于竞争劣势的物种被排除，从而促进物种共存(Tilman, 1994；Chesson, 2000)。此外，生境的保护是珍稀濒危植物保护最重要的环节之一(Volis, 2019；Caissy et al., 2020)，理解物种分布与环境因子的关联对于制定有效的保护措施至关重要(Austin, 2002; Pliscoff et al., 2014)。已有研究证明，非生物环境因子包括土壤养分、水分、地形等，以及生物环境因子，如物种多样性、邻体特征等，对目标物种的分布、生长、更新等过程均有不同程度的影响(Guisan & Thuiller, 2005；Elith & Graham, 2009；Chuine, 2010)。充分结合不同的环境因子，分析其与目标物种的关联，从而确定影响珍稀濒危物种种群动态的重要环境因子，可用于指导珍稀濒危植物的保护。

本研究通过在广东南岭国家级自然保护区建立20 hm永久监测样地，对样地内的植物进行调查监测，并测量土壤、地形等多种环境因子，分析该区域内华南五针松的种群现状、空间分布格局以及其生境特征。拟解决以下3个科学问题：(1)华南五针松个体的径级分布如何反映其种群发展趋势；(2)华南五针松个体具有怎样的空间分布格局；(3)哪些环境因子对华南五针松的个体数量分布具有重要的影响。

1 材料与方法

1.1 研究地点概况

广东南岭国家级自然保护区位于112°30′—113°04′ E、24°37′—24°57′ N，保护区面积为58 368.4 hm，位于广东省韶关市乳源县、清远市阳山县和连州市行政境界内，最高峰石坑崆，海拔1 902 m，为广东省第一高峰。该区域气候为亚热带温湿气候，年平均气温为17.7 ℃，降水量较充沛(3—8月为雨季)，年平均降雨量为1 705 mm，年平均相对湿度为84%。南岭保护区除保留了典型的亚热带常绿阔叶林外，还分布有沟谷雨林、针阔混交林、针叶林和山顶矮林等植被类型。南岭山脉在生物进化历史上具有重要的地位，既是古热带动植物的避难所，近代东亚温带、亚热带植物的发源地之一，也是现今我国16个生物多样性热点地区之一(班美玲等，2018)。

1.2 研究方法

1.2.1 20 hm森林动态监测样地建设 广东南岭国家级自然保护区20 hm样地于2020年开始建设，2021年5月完成样地调查。样地建设参照ForestGEO (https://forestgeo.si.edu)森林样地监测技术规范，样地建设主要包括样方设置、地形测绘、网格化设置和标记、以及树木标记。首先，用全站仪精确测定样地地形，确定样地范围并进行网格化设置，包括把样地分成500个投影面积为20 m × 20 m的基本样方单元，测量每个样方的地形，包括海拔(elevation)、坡向、坡度等，每个样方单元进一步分为16个5 m × 5 m的小样方。然后，对网格化样地进行界桩标记，包括对每个20 m × 20 m样方四角用花岗岩桩进行永久标记，对5 m × 5 m的小样方四角用PVC管进行标记；对样地内所有胸径≥1 cm的木本植物(含分枝)进行每木标记，标记内容包括画漆和挂牌。最后，进行每木调查，包括测量胸径、树高、冠幅，鉴定到物种，记录个体坐标等。

1.2.2 土壤理化性质测定和分析 参照ForestGEO森林样地土壤采样技术规范，选择固定点加随机点的方法确定采样点，以充分考虑土壤环境的空间异质性。把南岭大样地划分成238个30 m × 30 m的网格，取网格顶点为基准点；在每个基准点的东、东南、南、西南、西、西北、北和东北8个方向上随机选取1个方向，并在沿该方向2、4、12 m处随机取两个作为采样点，每个采样点安装PVC管并标记编号。总共设计700个土壤采样点，获取每个采样点0～15 cm的表层土壤，测定其水分含量(water，%)、pH值、有机质(OM，g·kg)、全钾(TK，g·kg)、速效钾(AK，mg·kg)、全磷(TP，g·kg)、速效磷(AP，mg·kg)、全氮(TN，g·kg)、硝态氮(AN1，mg·kg)和铵态氮(AN2，mg·kg)等土壤理化性质。土壤理化性质与华南五针松多度分布的关系分析以20 m × 20 m样方为单位。因此，采用普通克里格插值(ordinary krige)的方法分别对这10个土壤因子进行插值(Shen et al., 2019a)，得到样地范围内每个20 m × 20 m样方的土壤变量数据。

1.2.3 生物环境因子 以500个20 m × 20 m样方为单位，根据植物调查数据，计算每个样方的物种丰富度()、多样性指数(Shannon-Weiner index，′)，个体的密度(density)、平均胸径(DBH)和树高(height)，以及胸高断面积(basal area，BA)总和。

1.3 数据分析

1.3.1 20 hm样地华南五针松种群现状 基于南岭20 hm样地的调查和监测数据，对样地内的华南五针松种群现状进行分析，包括计算个体数量、个体大小、个体高度、径级分布、空间分布等基础信息。

1.3.2 空间点格局分析 基于个体的坐标，采用Ripley’s函数(Ripley’sfunction)分析华南五针松的空间分布格局(Ripley, 1977)。Ripley’s函数的定义为目标个体周围一定距离内的平均个体数量，通过函数计算华南五针松空间分布的观测值，并与随机模型产生的分布格局进行对比(de la Cruz et al., 2008)。Ripley’s函数通常转换为函数(function)。计算方法如下：

式中：()表示在某个距离内(= 100 m)华南五针松目标个体周围的同种个体的数量，与随机模型(模拟9 999次)生成的格局相比较；Π = π。

实际的function如果向上偏离包迹线，说明该物种是聚集分布，向下偏离包迹线则表示该物种为均匀分布，在包迹线范围内则为随机分布(Shen et al., 2013)。Ripley’s函数采用R语言(4.1.0)“spatstat”包来实现。

1.3.3 华南五针松多度与环境因子的零膨胀泊松模型 为了研究哪些生物环境因子和非生物环境因子影响华南五针松的分布与否(0到1)，以及分布数量的多少(1到N)，采用零膨胀泊松模型(Zero-inflated Poisson, ZIP)(缪柏其等，2008；赵静等，2020)。

以500个20 m × 20 m样方为单位，未出现华南五针松的样方为0，出现华南五针松的样方为(个体数量)。ZIP的拟合采用R语言(4.1.0)“pscl”包中的“zeroinfl”函数来执行。具体的公式如下：

zeroinfl(～ele+′+height+DBH++density+ba+water+pH+AN1+AN2+AP+AK+TN+TP+TK+OM|ele+′+height+DBH++water+pH+AN1+AN2+AP+AK+TN+TP+TK+OM，dist=“Poisson”)。

式中： 分隔符“|”将零膨胀模型分为两部分，分别对应二项过程和计数过程，每一部分使用的解释变量组合通过模型筛选来确定，获得最优模型。

具体方法如下：先建立饱和模型，即包含所有的17个预测变量，再通过不同预测变量之间的组合进行模型筛选，以赤池信息准则(akaike information criterion，AIC)来衡量模型拟合度，AIC越小表明模型拟合得越好，从中选出最优模型。模型筛选采用“MuMIn”包中的“dredge”和“model.sel”函数来完成(Shen et al., 2019b)。

为明确多重共线性的影响，最优模型采用方差膨胀因子分析(variance inflation factor analysis，VIF)来检验，自变量的VIF小于10则认为多重共线性较弱，适合建立回归模型(Ohlemuller et al., 2006)。

2 结果与分析

2.1 南岭20 hm2样地概况

广东南岭国家级自然保护区20 hm森林动态监测样地的海拔为1 008～1 254 m，是华南五针松种群的集中分布区域。样地的森林类型主要为亚热带常绿阔叶林，并伴有少量的落叶树种。样地共调查和监测了132 905个胸径≥ 1 cm的木本植物，属于63科127属的229个物种。优势物种多为樟科(Lauraceae)、壳斗科(Fagaceae)、蔷薇科(Rosaceae)、五列木科(Pentaphylacaceae)、山矾科(Symplocaceae)、杜鹃花科(Ericaceae)和冬青科(Aquifoliaceae)的物种。个体数量最多的物种为五列木()，多度为13 495；最大的个体为甜槠()，胸径达125.2 cm。

2.2 华南五针松径级及空间分布

华南五针松在样地内分布的多度为402，多数个体为大树，平均胸径达到17.5 cm，其中最大的个体胸径达到70.1 cm，径级分布近似钟型(bell-shaped)(图1)，中间径级的个体数量最多，大径级和小径级的个体数量都相对较少。所有个体的平均高度为7.6 m，高度范围为2～18 m。其在样地内主要的分布区域为西南角及东侧的海拔较高且地形比较陡峭的山坡和山脊上(图2)，尤其是在样地的西南角分布最为集中，在海拔较低且地势较平缓的样地中部(地形为沟谷)则几乎无分布，在空间上呈现出较为明显的聚集分布，以及较强的生境特异性(依赖于地形)。

图 1 南岭20 hm2样地华南五针松径级分布Fig. 1 Size-class distribution of Pinus kwangtungensis in the Nanling 20 hm2 forest dynamics plot

图 2 南岭20 hm2样地华南五针松空间分布Fig. 2 Spatial distribution of Pinus kwangtungensis in the Nanling 20 hm2 forest dynamics plot

采用函数对华南五针松的空间分布格局进一步的分析发现，其空间分布的观测值(黑线)远高于9 999次随机模型产生的包迹线(虚线)(图3)，并且随着距离的增大，实际的观测值与随机模型值的差异逐渐增大，图3结果显示，华南五针松在空间分布上为显著的聚集分布，即目标个体周围同种个体的数量远远高于随机模型的结果，且随着距离的增大这种聚集效应也在不断增强。

黑色实线为L function的实际观测值，灰色虚线为随机模型产生的包迹线的上限和下限。观测值向上偏离包迹线说明该物种是聚集分布； 向下偏离包迹线则表示该物种为均匀分布； 在包迹线范围内则为随机分布。The solid line is the observed L(r) and the dashed lines represent the simulations envelopes under the null models. Significant associations occur at scales where the observed L(r) exceeds the simulated envelopes. The observed positive departure of L(r) from the envelopes indicates an aggregated distribution, and a negative departure of L(r) from the envelopes indicates a regular distribution, within the envelopes indicates a random distribution.图 3 南岭20 hm2样地华南五针松空间点格局分析(L function)结果Fig. 3 Results of the spatial point pattern analysis (L function) for Pinus kwangtungensis in the Nanling 20 hm2 forest dynamics plot

2.3 华南五针松生境特征

采用零膨胀泊松回归模型可以研究生物和非生物环境因子对华南五针松分布的影响。零膨胀泊松回归模型分为两部分，包括二项式逻辑回归和泊松对数回归，前者解释的是环境因子如何影响华南五针松的有无(从0到1的过程)，后者解释的是华南五针松的多度分布(从1到N的过程)。通过不同变量之间的组合进行模型筛选，得到最优模型(表1)。最优模型的AIC最小，为782.2；权重(weight)最大，为0.349(权重从0到1，越大表示模型拟合度越好)。所有自变量的VIF值均小于10(1.53～7.40)，自变量之间多重共线性程度较低。

表 1 生物和非生物环境因子影响华南五针松分布的最优零膨胀泊松回归模型Table 1 Parameters of the optimal Zero-inflated Poisson regression model for the effects of environmental factors on the distribution of Pinus kwangtungensis

二项式逻辑回归的结果显示，土壤养分对于华南五针松的分布与否具有较密切的关联。与无华南五针松分布的区域相比，有华南五针松分布的区域通常具有较高的铵态氮(AN1)和速效钾(AK)，但其有机质含量(OM)却较低。海拔和生物环境因子在二项式逻辑回归中均不显著。从泊松对数回归的结果来看，海拔、土壤养分和生物环境因子均对华南五针松的多度分布有影响。在华南五针松分布的区域， 海拔越高其分布的个体越多；而土壤养分的结果则相反，华南五针松个体分布较多的区域，其土壤养分较低，包括铵态氮、总磷(TP)和总钾(TK)均与华南五针松的多度呈显著的负相关关系。从生物环境因子来看，华南五针松个体主要分布在较为成熟的林分中，如物种多样性(H’)较高且平均胸径(DBH)较大的区域，但其多度与周围植被的高度(height)却呈显著的负相关关系。

3 讨论与结论

气候环境的变化以及人为干扰的加剧，导致华南五针松种群数量急剧下降，目前原生的华南五针松种群已不多见，且分布的面积都较小，IUCN亦将其列为近危(NT)等级。因此，探索该物种的种群特征、分布以及影响其分布的环境因素和生态学过程，对华南五针松的保护和人工繁育等都具有现实的意义。本研究通过大规模的调查和监测，以及结合多维度的环境因子，探讨了华南五针松空间分布格局产生的生态学过程，为该物种的保护提供了重要的科学依据。

3.1 华南五针松“钟形”径级结构-种群衰退

种群径级分布结构是预测种群稳定性的重要指标。比如，个体数量随径级增大逐步减少，其径级结构表现为“倒J型”，通常认为种群是较稳定的。但是，幼树的减少将导致种群的径级结构呈现出“钟形”(叶万辉等，2008；Shen et al., 2013)，这种径级分布可能是种群更新困难和衰退的信号(Condit et al., 1998; Venter & Witkowski, 2010)。本研究结果显示，华南五针松的径级分布形状近似于钟形，小径级和大径级的个体数量均较少，意味着该种群可能存在更新缓慢的问题，一定程度上反映了种群衰退的迹象。此外，从南岭20 hm样地的600个1 m × 1 m 的幼苗监测样方的结果来看，2021年调查数据显示，林下无华南五针松幼苗的更新，进一步证实了该区域的华南五针松存在更新困难的现象，需要进一步加强监测和保护以防止种群的急剧衰退。

3.2 个体聚集分布种内竞争加剧

从空间分布来看，华南五针松的分布呈现较强的聚集分布，其偏好的地形为海拔较高且地形较为陡峭的山坡和山脊。通过空间点格局分析进一步证实，华南五针松的分布格局为强的聚集分布，且随着距离的增大聚集效应增强。物种聚集分布在热带和亚热带森林中很常见(Manabe et al., 2000; Mitsui & Kimura, 2000; Plotkin et al., 2000)，通过空间点格局的分析可以研究产生该格局的生态学过程(Janzen, 1970；Connell, 1971；He & Duncan, 2000)，是理解和模拟生物多样性空间变化的基础(He et al., 2002；Wright, 2002；Wills et al., 2006)。已有的研究认为，多种生态学过程，包括生态位隔离(Pielou 1961)、生境异质性(Harms et al., 2001)、邻体竞争(He et al., 2000)和扩散限制(Hubbell, 2011)等，均可以导致物种的聚集分布。

对于热带森林中的树种，现有的研究认为，扩散限制可以使不同物种在空间上形成隔离(种内聚集分布)，减少种间竞争排斥，如同种空间聚集的范围和程度与种子传播方式相关(Seidler & Plotkin, 2006)，通过动物传播种子的物种比通过风或喷射性传播的物种传播距离更远，通过动物传播的物种比风传播或喷射传播的物种聚集程度更低(Li et al., 2009)。松科物种的种子主要依靠风力来传播，仅有少量的种子通过鸟类或者啮齿类传播。因此，华南五针松种子传播的距离有限，一方面可以避免其与其他物种产生强烈的竞争(Tilman, 1994；Chesson, 2000)；另一方面导致其种内的竞争加剧，且种群的扩散受到限制，不利于种群的恢复和发展。

3.3 个体有无及数量多少依赖于不同的环境因子

生境异质性一直被认为是影响物种分布的主要因素之一(Li et al., 2009)。物种基于生态位分化的生境特异性是导致不同物种在不同生境中具有自身优势的原因(Harms et al., 2001)。物种聚集分布在适宜的生境中，使其具有更高的适合度(Cody, 1991)。比如，本研究关注的华南五针松，其适宜的生境主要是海拔较高、地形陡峭的山坡或者山脊(王献溥和李信贤，1989；陶翠等，2012)，而在地势较为平缓的沟谷中几乎没有分布，其分布格局表现出极强的生境特异性。

虽有研究报道华南五针松耐贫瘠(丁忠江等，2000)，但本研究结果却发现，与没有华南五针松分布的区域相比，华南五针松分布的区域具有较高的硝态氮和速效钾含量，说明其亦倾向于选择土壤养分较高的区域生长，以满足其基本的养分需求，这些土壤养分元素对于热带和亚热带森林树种的生长具有非常重要的影响(Wright et al., 2011；Santiago et al., 2012；Shen et al., 2014)。但是，其分布的区域具有较低的有机质，可能是由于其生境地形陡峭，枯枝落叶层更容易被风和降雨冲刷，导致腐殖质层流失(Garciaoliva et al., 1995; Shen et al., 2014)，土壤表层的有机质减少。针对华南五针松多度的分析(泊松对数回归)则得到更为不同的结果。首先，海拔的重要性显现，华南五针松在样地海拔范围内(1 008～ 1 254 m)均可以生长繁殖，海拔较高且地形陡峭的区域其多度要高于低海拔区域。然后，土壤养分的效应发生变化，高土壤养分的区域反而不利于多度的增加，这种现象可能与竞争的强度相关。比如，华南五针松为显著的聚集分布，种内竞争较强，且在环境较好的区域，资源的竞争加剧(Katabuchi et al., 2012; Shen et al., 2016)，竞争导致同种个体减少，群落结构趋异(Spasojevic & Suding, 2012)。最后，群落的结构(生物环境因素)对于华南五针松的多度分布产生了影响。本研究结果显示，华南五针松个体多分布在较为成熟(物种多样性高，个体较大)的群落中，作为阳生树种(沈燕等，2016)，冠层太高则不利于华南五针松个体的生长。华南五针松为群落的先锋和建群种之一，建群种在群落中处于优势地位(叶万辉等，2008；Shen et al., 2013)。因此，华南五针松同样在较为成熟的群落中处于优势的地位，而其多度与树高呈现出显著的负相关，可能是其生境多为陡峭的地形，该地形下群落冠层的高度受到限制。

本研究仅在局域尺度上研究了华南五针松的空间分布格局和生境特征，预测不同区域的格局时应考虑其他重要的环境因子的影响，如温度和光照等。此外，对于华南五针松的保护还需要结合实验的手段研究其生理、繁殖特征等，同时促进就地和迁地保护。