地果花粉供应格局成因和适生区分析

袁 川, 谭 露,2, 徐亮, 童 鑫,游章强, 陈 艳,*

1 绵阳师范学院生态安全与保护四川省重点实验室,生命科学与技术学院, 绵阳 621000 2 四川农业大学小麦研究所, 成都 611130 3 四川省环境监测总站, 成都 610041 4 华东师范大学生态与环境科学学院,浙江天童森林生态系统国家野外科学观测研究站,上海市城市化生态过程与生态恢复重点实验室,上海 200241

植物是否能适应环境在很大程度上取决于其有性繁殖对环境变化的响应速度[1-2],其中花粉供应起着关键作用[3-4]。花粉供应由性系统的表达和分配决定,众多案例揭示,被子植物的性别表达与分配具有环境依赖性[3,5- 7],在不同环境条件下被子植物会权衡对雌、雄的资源分配方式,使器官或功能在不同环境下更合理[8- 10]。在资源贫乏的环境中,植物往往会减少对雌性器官的投入,增加对雄性器官的投入[11]；如在高山环境中,随着海拔的增加,环境胁迫增加,植物对雄性功能的投入增加,而对雌性功能的投入减少[12-13]。环境条件的变化同样会影响动物的群落组成和行为,约有80%的被子植物依赖昆虫传粉[14-15],传粉者对植物性表达和性分配的影响已在许多植物类群中被证实[16- 18],植物可快速响应传粉昆虫种群动态和行为变化而形成不同的花粉供应模式和花部特征[19- 21],其中以依赖专一性传粉昆虫的植物对传粉者种群的变化尤为敏感[3]。

桑科榕属(Ficus, Moraceae)是被子植物中最大的属之一[22-23],其最典型的特征是具有隐头花序(即榕果),小花苞埋在隐头花序内部,花序顶部的苞片口是榕果与外界交流的唯一通道[24]。这种结构需要依靠专性传粉小蜂为其传粉,才能产生种子；而隐头花序也为传粉小蜂提供产卵、繁殖后代的场所[22-23],使榕-传粉榕小蜂(Agaonidae, Chalcidoidea)成为目前已知的专一性最强的互惠共生系统之一[25]。传粉榕小蜂有主动和被动两种不同的传粉行为。主动传粉小蜂具有花粉筐和花粉刷的特殊结构,有主动收集花粉和授粉的行为,具有很高的传粉效率；而被动传粉小蜂没有特殊的结构和主动传粉的行为,仅靠体表的毛发结构粘附花粉,花粉在小蜂随风扩散的过程中损失量较大,传粉效率较低。大量研究显示,植物花粉的相对供应量反映了物种的花粉传递效率,可采用花粉胚珠比(pollen-to-ovule ratio)来表征,在榕属物种中通常采用花药胚珠比(anther-to-ovule ratio,A/O)来衡量[18,20]。Kjellberg等[26]对142个榕属物种的传粉模式和花药胚珠进行了比较,发现两种传粉模式间有显著的分化,主动传粉模式的榕花药胚珠比多小于0.16,而被动传粉模式的榕花药胚珠比多大于0.23。

地果(FicustikouaBur.)是我国西南地区常见的榕属物种,其传粉小蜂隶属Ceratosolen属[25],具有主动传粉小蜂的典型结构和传粉行为[27],但Deng等[27]对41个自然种群的调查显示,地果的花药胚珠比有巨大的变幅(0.096—10.000)。同时,雄花分布格局存在差异,即有雄花集中分布在苞片口的OS类型(Ostiolar style)、雄花散生在果腔内的SS类型(Scattered style)及中间类型(IS型, Intermediate style),其中OS型的花药胚珠比显著低于SS型。不同类型呈现出明显的地域分布格局,四川盆地及周边地区主要为OS型,云贵高原主要为SS型,部分种群中存在多种类型共存的现象[27]。对大量榕属物种的研究显示,同一榕属物种的花药胚珠比很稳定[20,26],那么在地果中为何呈现出如此大的差异,环境条件在其中起到怎样的作用？地果为铺地藤本,耐受贫瘠环境,是良好的护坡植物,具有较好的生态修复应用前景。那么不同类型的地果适生区是否重叠？地果种群在不同适生区之间的扩散和定居是否受到环境条件的限制？

本研究依据地果分布样点数据,采用最大熵模型(MaxEnt)模拟不同类型地果的潜在分布区,提取主导环境因子,分析如下问题,旨在探讨环境条件对榕属物种花粉供应的影响、并为地果在不同区域的生态应用提供参考：(1)地果在中国的潜在分布区是否存在区域差异？(2)各分布区域的环境主导因子是什么？(3)各分布区域的环境条件与地果不同雄花分布类型和花粉供应水平之间有什么联系？

1 材料和方法

1.1 物种分布模型

物种分布模型(Species distribution models,SDMs)广泛应用于物种潜在分布区模拟及其对气候变化的响应、生物多样性保护等领域[28- 30]。其中,最大熵模型可通过部分信息预测物种分布概率,寻找关键环境贡献因子[31],对小样本具有较高的耐受性[32]。采用Phillips等[31]编写的MaxEnt程序,预测地果在中国的潜在适宜生存概率点,软件版本为MaxEnt 3.4.1。

1.2 数据来源

查询全球生物多样性信息数据库(GBIF, http：//www.gbif.org/)、中国数字植物标本馆(http://www.cvh.ac.cn/)、结合我们2013—2016年实地调查,得到514个地果自然分布点。中国植物志记载地果分布于我国西南、印度东北、越南北、老挝,但查询到的标本信息国外仅有8个点,其中夏威夷7个,因夏威夷不是地果的原生地,我们去除了8个国外分布点。采用系统重采样法修正取样偏差[33-34],即利用ArcGIS软件给分布点建立10 km缓冲区,每个缓冲区内随机保留一个位点,最终保留339个分布位点。

四大类31个环境因子数据均从公共数据库下载(表1)。数据分辨率统一转换成30″(约1 km×1 km),投影坐标系和地理坐标系分别设为PCS_WGS_1984_UTM_Zone_47N和GCS_WGS_1984,通过中国行政区范围(来自http://www.resdc.cn/)进行裁剪,并转换成ASC格式。

1.3 环境因子筛选与重要性评估

依据变量间的相关性和变量对模型的贡献筛选环境因子,保留相关性大于0.85的变量中贡献更大者。环境因子之间的相关性采用ENMTools软件(http://purl.oclc.org/enmtools)计算[35]。

环境因子对MaxEnt模型的相对贡献常采用刀切法(Jackknife)结果、贡献率(Percent contribution)和置换重要性(Permutation importance)三个侧重点不同的指标进行评价[36]；分别赋予贡献率和置换重要性各50%的权重,计算变量的综合得分。三个指标在MaxEnt程序默认函数模型,10倍交叉验证(10-fold cross-validation),模型的最大背景点数量为10000,最大迭代3000次条件下预测得到。

1.4 预测模式筛选

使用MaxEnt程序中的5种函数模式(线性函数L、二次函数Q、交互函数P、常数分段函数T、线性分段函数H)及其组合共计13种模式(L、Q、H、T、LQ、HQ、LQP、LQT、QHP、QHT、QHPT、LQHPT、Auto)筛选最优函数模式；使用样本数据中随机的75%作为训练数据,剩余25%作为测试数据,3个重复进行模拟。基于小样本赤池信息量准则(AICc)对13种模式在调控系数分别为0.5、1、2和4的情况下进行测试[37]。使用ENMTools[35]中的方法计算AICc值,以AICc值最低的模型为预测模型。

1.5 模型预测与评估

将分布点和环境因子数据导入MaxEnt 3.4.1构建当前潜在分布模型,使用受试者工作特征曲线(Receiver operating characteristic,ROC)评估模拟结果的精度,通过曲线下面积(Area under the curve,AUC)表示,AUC值越接近于1,模型表现效果越好,使用刀切法计算单个环境因子重要性,使用10倍交叉验证法评估模型的有效性。

根据筛选结果选择函数模式,设定模型调控系数再次进行模拟,计算变量的贡献。根据Liu等[38]提出的阈值选择原则,结合Radosavljevic等[39]和Anderson等[40]的研究结果,选择训练数据的10%存在率(10 percentile training presence)作为大样本情况下的阈值,建模数据敏感度与特异性之和最大的逻辑阈值(Maximum training sensitivity plus specificity logistic threshold)作为小样本情况下的阈值。使用ArcGIS结合阈值进行分类,得到地果的适宜分布区。

表1 环境因子及数据来源

采用相同的方法对地果的总体和区域分布分别进行模拟。总体预测以修正后的339个分布点为依据,将预测结果图划分为10 km×10 km的栅格,提取适宜生境的占有比例。获得占有率高于90%的连续分布区域2块(四川盆地,贵州高原),占有率低于50%的破碎化区域1块(横断山区-云南高原),分别以3个区域中地果分布点为依据,进行区域模拟。

1.6 限制因子分析

限制因子是致使模型预测值产生最大变化的变量,采用多元环境相似度面分析(Multivariate environmental similarity surface, MESS)确定驱动MESS值的变量、最不相似变量分析(Most dissimilar variable, MOD)寻找最不相似变量,进行限制因子的提取和映射。所有分析使用maxent.jar中的density.Novel工具完成[41]。MESS负值是外推的结果,表示至少有一个变量的值超出了参考集的范围,用红色表示；正值表示内推结果,值的大小表示正常程度,100表示完全正常,用蓝色表示。一般情况下,模型的内推结果比外推结果可信度更高[41]。

1.7 区域间环境因子的评价与比较

综合贡献率与置换重要性得分,提取主要环境因子及其与地果适宜性的相关性或地果生存最适值信息；分别提取四川盆地 60个、横断山区-云南高原72个、贵州高原 77个自然分布点的关键环境因子,利用Kruskal-Wallis检验判断各区域环境因子的差异。

1.8 地果各类型的分布与环境因子

Deng等[27,42]按照地果雄花分布格局,将地果区分为不同的类型,将Deng等[27,42]的种群投射到MaxEnt的模拟图中,共56个种群(表2),其中OS种群仅有OS类型、SS种群仅有SS类型、Mix种群具有多种类型的榕果,计算种群中OS、SS类型出现的相对频率及花药胚珠比。榕果的花药胚珠比采用Deng等[27]的数据,计算方法为,A/O=雄花数/雌花数×2(每朵雄花有两个花药)[26]。不同类型地果的花药胚珠比差异采用Kruskal-Wallis检验。

采用卡方检验(Chi-square test)分析地果三类种群在三个区域是否随机分布,根据分布区域中OS、SS、Mix种群的出现频度生成3×3二联表进行检验,重复2000次获得p值。

提取各种群的关键环境因子,使用广义线性模型(Generalized linear model, GLM)检验三类种群环境因子的差异,利用斯皮尔曼相关性检验(Spearman correlation test)检验OS、SS类型出现频率与环境因子的相关性。

所有检验都在R3.6.2软件中进行。

表2 各区域取样种群及其地果类型

区域Regions种群Populations经度Longitude/(°)纬度Latitude/(°)海拔Altitude/m类型Type鸣鹫103.75 23.42 1979.00 SS普者黑104.05 24.09 1538.00 SS师宗104.18 24.78 1760.29 SS武定102.39 25.54 1795.00 SS砚山104.29 23.65 1611.00 SS曲靖103.71 25.54 1900.00 Mix宣威104.06 26.22 2033.00 Mix

2 结果

2.1 模型表现力评价与适生区分布格局

总体预测使用15个环境因子、LQT的组合模型、调控系数2；10倍交叉验证训练集(Training data)和测试集(Test data)的平均AUC值都大于为0.945,说明预测结果的可信度很高。

3个分布区域的优选预测条件及可信度如下：四川盆地14个环境因子、函数模式Q、调控系数0.5,训练集和测试集的平均AUC值为0.994；贵州高原14个环境因子、函数模式LQ、调控系数0.5,模拟训练集和测试集的平均AUC值为0.987；横断山区-云南高原14个环境因子、函数模式Q、调控系数2,预测结果的平均AUC值为0.979。

总体预测根据训练数据的10%存在率切点0.3372划分潜在适宜和不适宜分布区,地果集中分布在我国西南地区(图1a)。区域预测使用建模数据敏感度与特异性之和最大的逻辑阈值作为划分潜在适宜和不适宜区的切点,区域预测的潜在适应分布区保持异域分布,说明三个区域的环境条件有较大差异(图1b)。

图1 地果适宜分布区模拟图Fig.1 Simulated suitable distribution area of Ficus tikoua图a为总体预测结果及划定的三个区域生境；图b为区域预测结果,不同颜色代表不同区域的预测结果；点代表地果种群的类型

2.2 区域分布的环境贡献因子

关键环境因子的贡献率、置换重要性、综合得分及最适量见表3。四川盆地种群的适合度与海拔、温度季节变化标准差、昼夜温差月均值等环境因子呈负相关(图2)。昼夜温差月均值、年潜在蒸发量、最干季度降水量在贵州高原种群中呈现出单峰适应性曲线(图2)。横断山区-云南高原种群的适合度与温度季节变化标准差呈负相关,而与太阳辐射量呈正相关,气候湿度指数在种群中呈单峰适应性曲线(图2)。

图2 不同区域地果种群的主要环境因子适宜性Fig.2 The suitability of key environmental factors in the different distribution regions of Ficus tikouaPET: 潜在蒸发量Potential evapotransporiation

2.3 区域间扩散的环境限制因子

分别对三个区域的预测结果进行多元环境相似度面分析,结果显示地果潜在适宜区都分布在正值范围内(图3a、c、e),均为内推结果,模型精度较高。

利用最不相似变量分析提取各区域地果种群向区域外扩散的限制因子(图3 b、d、f)。海拔、降水量及其变化限制了四川盆地种群向高原的扩散(图3b、图4)；降水差异、太阳辐射量限制了贵州高原种群向外扩散(图3d、图4)；等温性、温度的季节变化限制了横断山区-云南高原种群向其他区域的扩张(图3f,图4)。

表3 影响地果在各模拟区域分布的主要环境因子及其贡献

2.4 区域间关键环境因子的差异分析

三个区域主要贡献环境因子及区域间扩散的限制性因子总结如图4、这些环境因子在三个区域的差异如表4,重要因子在区域间的差异显著,四川盆地海拔低、温度的季节变化大、昼夜温差小；横断山区-云南高原海拔高、太阳辐射强、温度季节变化小、昼夜温差小；而贵州高原气候湿度指数高、降水变异系数小、最冷季与最干季降水量高。

表4 区域间的关键环境因子差异

图4 区域分布的贡献因子与区域间扩散的限制因子Fig.4 Key environmental factors contributing to the regional distribution or restricting expanding among regions圈内表示主要贡献因子,圈外为扩散限制因子,虚线代表扩散方向

2.5 地果不同类型的区域分布及其环境因子差异

地果的雄花分布格局在不同地域呈现出显著差异,卡方检验显示OS、SS、Mix三类种群非随机分布(x2= 22.43,P<0.001),其中OS种群主要分布在四川盆地,Mix种群主要分布在贵州高原,而SS种群主要分布在横断山区-云南高原(图1b)。SS类型的花药胚珠比显著高于OS类型(图5)。

提取三种类型所有种群的环境因子,海拔、温度季节变化标准差、最冷季降水量、太阳辐射、最冷月最低温在三类种群的分布有显著差异(表5),并与四川盆地、横断山区-云南高原、贵州高原的环境因子差异表现出一致的格局。OS种群分布于低海拔、温度季节变化明显、最冷季降水量较小、最冷月最低温较低的四川盆地区域；SS种群分布于高海拔、强太阳辐射、温度季节变化较小的横断山区-云南高原区域,而Mix种群分布于环境条件居中的贵州高原区域。检验各自然种群中OS、SS类型的出现频率与环境因子的相关性也证实,OS型的分布与海拔显著负相关(Rho=-0.469,P<0.001)、与温度季节变化标准差显著正相关(Rho=0.529,P<0.001)、与最冷季降水量显著负相关(Rho=-0.296,P=0.029)；SS型的分布与海拔显著正相关(Rho=0.494,P<0.001)、与温度季节变化标准差显著负相关(Rho=-0.527,P<0.001)、与最冷季降水量显著正相关(Rho=0.283,P=0.038),但与太阳辐射量(Rho=0.212,P=0.123)、最冷月最低温(Rho=0.220,P=0.120)没有显著的相关性。

3 讨论

3.1 环境条件对不同类型地果分布格局的影响

地果的生长环境在不同地域呈现出不同的特征,环境条件的差异限制了地果种群在区域间扩散。MaxEnt模拟显示,地果非连续分布,可分为三个生境区域(图1a)。比较区域的主要环境贡献因子和扩散限制因子、以及地果实际分布点的主要环境因子,结果均显示三个区域的环境条件有明显的分异。四川盆地与横断山区-云南高原区域呈现出相反的环境特征,四川盆地海拔低、温度季节变化大,横断山区-云南高原海拔高、温度季节变化小、太阳辐射强；而贵州高原的环境条件介于两者之间(表3、表4)；三个区域的环境条件差异限制了区域间的扩散(图3,4)。

地果不同类型的分布呈现出环境依赖性。将Deng等[27,42]的56个自然种群的类型投射到三个生境区域中,卡方检验显示不同类型并非随机分布,OS类型主要分布于四川盆地,SS类型主要分布于横断山区-云南高原,而贵州高原则主要为多种类型的共存(图1b)。

表5 三类种群的环境因子差异

图5 苞片口型(OS)和果腔散布型(SS)地果的花药胚珠比 Fig.5 Anther-to-ovule ratio of ostiolar style (OS) and scattered style (SS) Ficus tikoua types

地果不同类型的自然种群呈现出与模型模拟预测一致的区域环境条件差异。尽管OS、SS种群(仅有一种类型的种群)的分布并不仅局限在一个生境区域,比较56个自然种群的环境因子仍揭示出与预测区域一致的环境差异,OS种群分布的海拔、太阳辐射显著低于SS种群、而温度季节变化标准差显著高于SS种群(表5)；检验OS、SS型果在种群中的出现频率与环境因子的相关性也证实,OS型的出现频率与海拔呈显著负相关、与温度季节变化标准差呈显著正相关,而SS型正好相反；进一步揭示地果不同类型可能的环境适应性。

3.2 环境条件对地果性别分配策略的影响

被子植物的性表达和性分配具有环境依赖性,往往呈现出适应环境的可塑性[6,43],主动传粉的榕花药胚珠比显著低于被动传粉榕,且同一榕属物种的花药胚珠比极为稳定[26],而地果的自然种群中呈现出花药胚珠比的巨大变幅,其SS型的花药胚珠比显著高于OS型(图5)。比较不同类型地果的适生区发现,其环境条件显著不同,OS型分布于低海拔、温度季节变化大的区域,而SS型分布于高海拔、温度季节变化小、太阳辐射强的区域,即在高海拔、温度季节性变化小、太阳辐射强的区域地果供应给传粉小蜂的花粉量显著增高。

环境因子可能通过影响地果传粉小蜂,进而对地果的花粉供应格局产生选择效应。海拔是一个综合的环境因子,在高山环境中,海拔的增高往往伴随着温度的降低、生长季节的缩短和有效积温的降低；这些环境胁迫可直接作用于植物,也可通过影响传粉者间接影响植物。许多植物响应环境胁迫的增加和传粉者的稀少,随着海拔的升高而增加对繁殖器官和雄性功能(花粉)的投入[12,44],但地果的OS、SS类型分布于不同纬度,分布区海拔的增高并未伴随温度的降低和生长季的缩短。比较地果OS、SS种群分布区的年均降水、最冷月最低温显示(表5),分布于高海拔的SS种群因纬度更低,其水热条件都不低于OS种群,因此海拔对地果性表达的影响可能是通过传粉小蜂完成。

高海拔、强太阳辐射可能降低地果传粉小蜂成体的寿命和活力,降低了花粉的传递效率,迫使宿主地果增加对花粉的投入。已有研究证实,海拔对传粉者有较大的影响,随着海拔的增高传粉者的多样性、丰富度和活动能力都会降低[45-46]。地果为主动传粉的榕属物种,Deng等[27]对38个自然种群的地果传粉小蜂的形态扫描显示,小蜂的关键结构在不同分布区域并没有明显变化,都呈现出主动传粉榕小蜂的典型形态特征。主动传粉具有很高的传粉效率,有利于宿主植物的繁殖,但对于传粉小蜂而言,却需要消耗大量的能量,可能会影响雌蜂对产卵的能量投入。传粉榕小蜂成体体型小(成体3—5 mm),出飞后就不再进食,一般只能存活1—2天[47],因此传粉小蜂成体的活力和寿命可能会影响雌蜂对产卵和授粉行为的权衡,进而影响雌蜂的传粉效率。太阳辐射特别是紫外辐射会降低榕小蜂的活力和寿命[48],SS种群的太阳辐射显著高于OS种群,这可能会降低传粉小蜂成体的活力,迫使雌蜂将更多的能量投入到产卵过程,而降低传粉投入。适应于传粉者传粉投入的降低,地果则增加对花粉的投入。Deng等[27]对地果传粉小蜂体表的扫描显示,在高海拔的云南高原地区,地果传粉小蜂体表携带的花粉要高于低海拔的四川盆地地区,这是否意味着地果传粉小蜂的传粉行为在不同区域有细微的变化还需要进一步的实验加以揭示。模型模拟同时显示,横断山区-云南高原区域的适宜性随太阳辐射的增强而增加,这可能是因为太阳辐射增强能增加温度,温度升高会增加化合物的挥发性,榕属物种依赖于特定的挥发性化合物吸引传粉小蜂,化合物挥发性的增强有利于吸引传粉者[49-50]。

温度的季节变化可能通过影响地果繁殖物候的株间同步性、地果传粉小蜂种群的季节动态,而影响地果的性分配方式。榕与传粉榕小蜂是目前已知的专一性最强的共生体系之一,因传粉榕小蜂成体成活时间很短,若在传粉小蜂出飞期缺乏接受期的榕果将导致共生关系的解链,因此株间不同步、全年产生榕果的繁殖物候在榕属物种中较为常见[22]。地果OS类型分布区的季节性强于SS类型分布区,前者的温度季节变化标准差显著高于后者,而等温性的差异限制了后者向前者区域的扩散(图4)。被子植物的物候对温度变化较为敏感[51],温度的季节性往往会增加繁殖物候的同步性[52],亚热带北缘榕属物种繁殖物候的株间同步性增强,秋冬季节的低温减少了新生榕果的数量,同时冬季的极端低温还会导致大量落果[53-54],使传粉小蜂种群遭受严重的瓶颈效应(榕果是小蜂幼虫生存和发育的唯一场所),种群急剧降低、甚至不能维持[55]。Zhao等[53]对四川盆地地区地果种群繁殖物候的观察显示,只有极少数的榕果能越冬；每年春季地果大量萌生新一代的榕果,以扩繁冬季急剧减小的传粉小蜂种群,但因越冬榕果极少,春季出飞的传粉小蜂数量较少,新一代榕果能获得传粉小蜂雌蜂的不确定性较高,因而每年春季有大量的新生榕果因没有授粉而脱落。冬季传粉小蜂种群的急剧下降、春季获得传粉小蜂的不确定性增高可能促使OS类型调整繁殖投入策略,减少对单只传粉小蜂的花粉供应,而持续产生新榕果以等待传粉小蜂雌蜂的进入。我们在绵阳的野外观察显示,地果在整个春季均持续产生榕果。比较OS、SS类型分布点的最冷月最低温显示(表5),OS类型分布区的冬季最低温显著低于SS分布区,但这是否导致两种类型的繁殖物候、传粉小蜂种群季节性波动等方面的差异,仍然需要进一步的实验和观测加以证实。

4 结论

使用最大熵模型,基于339个地果分布位点模拟地果的潜在适宜区,对不同分布区域的环境贡献因子、区域间扩散的环境限制因子进行分析,探究地果雄花分布类型和花粉供应的地域分布格局成因。

结果显示,地果的不同类型对应着不同的适生区,OS型分布于低海拔、温度季节变化明显的区域,而SS型分布于高海拔、强太阳辐射的区域。

SS分布区的高海拔、高太阳辐射,可能降低传粉小蜂雌蜂的活力和寿命,雌蜂权衡产卵和传粉,降低主动传粉的能量投入,SS类型可能通过增加对单只小蜂的花粉供应来补偿传粉小蜂传粉投入的降低。

OS分布区的强季节性增加了地果繁殖物候的株间同步性,冬季传粉小蜂种群急剧降低、春季可获得传粉小蜂的不确定性增加,OS种群减少对单只小蜂的花粉供应。

致谢：杭州师范大学董鸣教授、华东师范大学陈小勇教授帮助写作,邓竣尹帮助数据分析,王勤、付荣华帮助样点整理，特此致谢。