近红外–微区X射线荧光多模态表型分析揭示南美假山毛榉属全部物种的互补化学与元素变异

《Plants》:Multimodal NIR–Micro-XRF Phenotyping Captures Complementary Chemical and Elemental Variation Across All South American Nothofagus

【字体: 大 中 小 】 时间:2026年11月03日 来源:Plants 4.1

编辑推荐:

  对保守林木种质的快速表征可增强迁地保护,以应对全球变化对森林生物多样性日益加剧的压力。研究人员评估了近红外–微区X射线荧光(NIR–micro-XRF)表型分析能否在保存于同质园中的所有十个南美假山毛榉属(Nothofagus)物种(N. alessandri

  
对保守林木种质的快速表征可增强迁地保护,以应对全球变化对森林生物多样性日益加剧的压力。研究人员评估了近红外–微区X射线荧光(NIR–micro-XRF)表型分析能否在保存于同质园中的所有十个南美假山毛榉属(Nothofagus)物种(N. alessandrii、N. alpina、N. antarctica、N. betuloides、N. dombeyi、N. glauca、N. macrocarpa、N. nitida、N. obliqua和N. pumilio)及天然杂交种N. × leoni中解析互补的特征信号。研究人员结合了表型分析、元素分辨推断、重复交叉验证偏最小二乘判别分析(PLS-DA)、多模态融合以及无监督双向正交偏最小二乘(O2PLS)。预测分析聚焦于由各12株独立树木代表的九个物种(n = 108);N. nitida仅由一株树的三次技术测量代表,保留用于描述性分析但排除于预测建模之外。全光谱X射线荧光(XRF)优于近红外(NIR)(准确率0.822 ± 0.056对0.701 ± 0.064),而区块均衡融合获得最高准确率(0.842 ± 0.040),但较单独XRF仅略有提升。O2PLS显示45.8%的NIR变异和52.1%的XRF变异属于共享结构,其类群相关组织形式(PERMANOVA R2 = 0.605)强于传感器特异性组分。Fe、Zn、Ca和Mn相关信号显示出最强的类群分化。若能在独立收集和自然群体中得到验证,该多模态框架可支持活体收集的分类学验证、表型表征和保护审计。
**研究背景与意义**

活体收集是植物迁地保护的重要组成部分,但其保护价值取决于可靠的分类学身份及对所包含生物多样性的有效表征。这一多样性涵盖分类学、功能与进化等互补维度,而种内表型变异构成常被比较评估所低估的重要层面。对于长寿树木而言,来源信息不完整和代表性受限可制约收集评估与管理。高维光谱表型可通过提供分类学身份与表型变异的互补信息来应对这些挑战,前提是其生物学解释与预测性能保持区分。在此背景下,对分类学上全面的迁地收集进行可扩展的表型表征,可为保护面临快速变化环境的林木谱系提供基线。

叶片反射光谱构成整合性表型而非单一预设性状的测量。在近红外(NIR)和短波红外区域,变异源于涉及O-H、C-H和N-H键的宽泛合频与倍频吸收,以及与组织和颗粒结构相关的物理散射。因此,NIR光谱整合了水相关属性、碳水化合物、木质素、纤维素、蛋白质、脂质和氮相关化学的协变。叶片光谱可区分近缘类群、在种子植物中保留系统发育信号并在对比环境条件下恢复遗传身份。然而,光谱表型并非环境不变:其整合了基因型、环境及其互作。

X射线荧光(XRF)提供物理和生物学上截然不同的光谱表型。XRF不测量分子吸收和散射,而是检测内层激发后发射的元素特征次级X射线。植物离子组是由土壤和根际条件、选择性吸收、运输、储存及在遗传和发育控制下的稳态调节所塑造的整合表型。台式与便携式XRF系统可快速分析植物材料中的多种元素,但信号幅度取决于样品呈现、颗粒大小、厚度、几何形状和校准;因此标准化干燥、研磨和采集至关重要。XRF离子组学已揭示植物和标本材料中协调的元素模式,但全XRF光谱编码森林树木类群身份并补充NIR衍生分子组成和散射相关表型的能力仍相对未被充分探索。

以往植物研究主要将NIR和XRF结合用于定量叶面分析而非类群判别。Comino等结合了干燥研磨橄榄叶的FT-NIR和便携式能量色散XRF测量,评估了低层和中层数据融合用于养分预测和作物营养诊断。该工作确立了两模态的分析互补性,但未解决从同一标准化材料获得的配对NIR和微区XRF测量是否恢复共享与模态特异性类群相关结构的问题。据研究人员所知,该问题在长期迁地树木收集中仍未被充分探索。

长期建立的同质园为评估这种互补性提供了严格环境。采样自然群体将把持久的类群分化与对局部环境的即时可塑性响应相混淆。在同一位点的长期生长减弱但不消除当代环境异质性,并允许持久差异被解释为在共同条件下维持的表型分化而非纯遗传效应。南美假山毛榉属是保护相关测试系统,因其类群涵盖落叶和常绿策略、分布于对比鲜明的温带环境,并包括杂交、渐渗和质体捕获可模糊物种边界的网状进化历史。研究集包括受威胁的智利特有種Nothofagus alessandrii、N. glauca和N. macrocarpa。N. alessandrii尤其分布范围受限,存续于遭受反复火灾和Pinus radiata入侵的破碎化森林残存中。长期建立的活体收集因此对保护研究有价值,尽管来源记录不完整和种内代表性受限可限制其对野外多样性的捕获程度。

**主要关键技术方法**

研究人员采集了同质园中十个南美假山毛榉属物种及天然杂交种N. × leoni的同一匀质化叶片粉末,分别获取NIR反射光谱(NIRQuest512光谱仪,895.915–2514.925 nm,512波长变量,经Savitzky–Golay一阶导数+SNV预处理保留502预测变量)和微区XRF光谱(Bruker M4 TORNADO Plus,50 kV/600 μA,1–15 keV保留1400通道,log(count+1)变换)。采用重复分层三折交叉验证(重复20次,60个评估分区)的PLS-DA进行分类,并进行特征缩放早期融合与区块均衡融合。无监督O2PLS分解共享与传感器特异性变异。统计推断使用Welch ANOVA、Kruskal–Wallis检验、PERMANOVA、PERMDISP及pvclust多尺度自助法。样本队列来源为智利洛斯拉戈斯地区Puerto Octay的Parque VerdeNativo活体收集,树木已在同一位点生长约18–25年,叶片于2025年1月8–9日单次采集。

**研究结果**

2.1 NIR和XRF揭示截然不同、模态特异性的分类学结构

尽管对齐分析矩阵包含117个样本单元,但这些不代表117个独立生物学观察。预测分析使用主要平衡数据集中的108株独立树木(九个非杂交物种,各12株)和包含N. × leoni的探索性数据集中的114株独立树木。三个N. nitida光谱源自一株生物学树木,仅保留用于描述性类群水平分析。XRF和NIR产生物理上截然不同的光谱表型。XRF光谱以窄的元素相关峰为特征,最显著者低于9 keV,叠加于宽泛连续谱之上。相反,原始NIR光谱显示宽泛且强烈重叠的吸收特征以及个体间基线反射率和光谱幅度的显著变异。物种均值热图显示两模态均存在结构化的种间变异,但分类学排列不同。在NIR树状图中,N. × leoni与N. glauca聚类(AU = 100,BP = 100),N. alpina与N. pumilio聚类(AU = 100,BP = 100),N. antarctica与N. dombeyi聚类(AU = 100,BP = 100),N. macrocarpa与N. obliqua聚类(AU = 100,BP = 99)。在XRF树状图中,N. × leoni转而与N. obliqua聚类(AU = 99,BP = 99),而N. glauca-N. macrocarpa、N. alessandrii-N. antarctica和N. betuloides-N. dombeyi配对均获得AU = 100和BP = 100。若干较深的NIR节点支持度较低(AU/BP = 78/74),而两个主要的较深XRF分区得到强支持(AU/BP = 99/99)。因此,两模态均解析出强支持的终端亲缘关系,但这些亲缘关系和更高阶的表型结构是模态特异性的。在类群均值水平,叶习性解释NIR的18.8%、XRF的33.8%和特征缩放融合差异的29.8%。精确随机化支持XRF(R2 = 0.338,p = 0.0121)和融合(R2 = 0.298,p = 0.0121)的叶习性相关组织,而相应NIR模式较弱且未达0.05阈值(R2 = 0.188,p = 0.0545)。PERMDISP显示NIR(p = 0.412)、XRF(p = 0.406)或融合(p = 0.848)均无多变量离散度不等的证据。在物种均值水平,亚属归属(Lophozonia对Nothofagus,排除单型Fuscospora类群N. alessandrii)分别解释NIR、XRF和特征缩放融合差异的28.4%、38.9%和36.1%(p = 0.0156、0.0072和0.0080)。由于叶习性和亚属在类群内不变,这些结果描述的是宽泛的组水平光谱组织而非独立于分类学身份的效果。

2.2 元素分辨XRF变异为类群判别提供可解释基础

经分析质量控制和SNIP背景校正后,保留八个元素相关信号用于比较分析:Si、P、S、K、Ca、Mn、Fe和Zn。Na、Mg、Al和Cu的候选峰未充分分辨而被排除。Cl亦被排除,因其特征K壳X射线发射与X射线管的铑L壳发射重叠,妨碍可靠归属。Welch ANOVA表明经Benjamini–Hochberg调整后所有八个信号均存在类群效应。效应量最大者为Fe(η2 = 0.720)、Zn(η2 = 0.691)、Ca(η2 = 0.575)和Mn(η2 = 0.533),其次为Si(η2 = 0.487)、K(η2 = 0.466)、S(η2 = 0.384)和P(η2 = 0.198)。P的证据最弱(BH调整p = 0.0439),其类群效应在平衡九物种Welch敏感性分析中不显著(BH调整p = 0.0581)。相反,Si、S、K、Ca、Mn、Fe和Zn的类群效应得到Welch ANOVA、Kruskal–Wallis检验和平衡敏感性分析的一致支持。基于八个保留信号的PLS-DA模型选择六个潜成分,获得重复交叉验证准确率0.598 ± 0.058(κ = 0.551 ± 0.065)。在平衡九物种敏感性分析中,准确率为0.628 ± 0.057(κ = 0.581 ± 0.064)。Mn、Fe、K和Si为四个最具影响力的信号(相对重要性 = 100.0、97.1、88.5和84.7)。在60个评估分区中,平均类群敏感性最高者为N. antarctica(0.988),其次为N. obliqua(0.846)和N. alpina(0.817),但对N. betuloides(0.104)、N. glauca(0.067)和N. × leoni(0.000)极低。因此,缩减信号集保留了实质性但不完全的判别信息:在同一平衡九物种框架内,其准确率为0.628,而全XRF光谱为0.822。八信号模型相对于全光谱XRF的较低性能表明,保留的整合峰面积未捕获完整光谱中存在的所有类群判别信息。

2.3 全光谱XRF优于NIR,多模态融合带来增量增益

在主要平衡分析中,14成分全光谱XRF模型获得准确率0.822 ± 0.056(κ = 0.799 ± 0.063),而13成分NIR模型为0.701 ± 0.064(κ = 0.664 ± 0.072)。特征缩放早期融合选择21个成分,获得准确率0.833 ± 0.048(κ = 0.812 ± 0.054),而区块均衡融合选择31个成分,获得0.842 ± 0.040(κ = 0.822 ± 0.045)。共识准确率遵循相同排序:NIR为0.759,XRF为0.824,早期融合为0.870,区块均衡融合为0.889。分区水平对比显示相对于XRF的增益温和且异质。相对于XRF,平均准确率差异为早期融合+1.16个百分点、区块均衡融合+1.99个百分点。相应经验2.5–97.5百分位范围包括零(分别为?7.0至11.1和?7.0至13.9个百分点)。早期融合在45.0%的分区中超过XRF并在23.3%中持平,而区块均衡融合在53.3%中超过XRF并在13.3%中持平。由于评估样本在重复折叠间重叠,这些配对分布被描述性解释而非正式推断。类群水平敏感性进一步表明融合效果并不均一。相对于XRF,区块均衡融合提高了N. glauca(0.68至0.89)、N. macrocarpa(0.65至0.74)和N. obliqua(0.93至1.00)的平均敏感性,但降低了N. dombeyi(0.77至0.59)的敏感性。N. betuloides在各模态中始终是最难准确解析的物种(敏感性 = 0.50–0.56)。包含杂交种N. × leoni的探索性10类群分析保持与主要分析相同的性能排序,准确率分别为NIR 0.661 ± 0.065、XRF 0.775 ± 0.045、早期融合0.800 ± 0.045和区块均衡融合0.805 ± 0.054;相应κ值分别为0.622 ± 0.073、0.749 ± 0.050、0.777 ± 0.050和0.783 ± 0.060。平衡准确率显示相同排序(分别为0.631 ± 0.064、0.739 ± 0.045、0.773 ± 0.050和0.779 ± 0.058)。因此,纳入代表性不足的天然杂交种降低了总体性能但未改变四种建模策略的相对排序。

2.4 无监督O2PLS揭示大量共享NIR-XRF类群相关结构

为确定多模态互补性反映的是共同潜在表型而非仅监督分类性能差异,使用无监督O2PLS分解平衡九物种NIR和XRF区块。交叉验证选择四个联合成分、四个NIR特异性成分和三个XRF特异性成分。联合变异分别占标准化NIR变异的45.8%和标准化XRF变异的52.1%,而传感器特异性成分分别占22.1%和22.8%;其余部分为残差。配对NIR和XRF联合得分在四个成分间一致相关(Pearson r = 0.712、0.549、0.634和0.661),证明存在大量跨模态协变。尽管物种身份未用于估计O2PLS成分,类群相关组织在整合共享空间中最强(PERMANOVA R2 = 0.605,伪F = 18.96,p = 0.0001),相较于NIR特异性(R2 = 0.212,伪F = 3.33,p = 0.0001)和XRF特异性(R2 = 0.300,伪F = 5.32,p = 0.0001)空间。平均轮廓宽度在共享空间中适度为正(0.056),但在NIR特异性(?0.103)和XRF特异性(?0.099)空间中为负,表明存在结构化的类群位移而非完全分离为离散簇。尽管个体重叠,物种质心仍占据共享得分图的分化区域。复合共享载荷分布于多个NIR和XRF区域,支持共同的跨模态光谱表型以及额外的模态特异性结构。

**讨论与结论**

从保护角度看,核心发现是长期建立的活体收集保留了跨互补光谱域可恢复的连贯类群相关表型结构。从同一匀质化叶片材料获得的NIR和XRF测量解析出不同的类群水平排列和分类误差模式,表明两模态捕获了叶片变异的互补维度。全光谱XRF提供最强的独立判别,得到多个元素相关信号实质性种间差异的支持。NIR对分类学分配准确性较低但贡献了非冗余的分子组成和散射相关信息,而低层融合产生最高的总体点估计。然而,相对于单独XRF的增益温和、在验证分区间异质且依赖类群。更广泛的启示因此不是增加预测变量维度必然改善分类,而是物理上截然不同的光谱域恢复表型分化的互补组分,其有用性取决于类群和分析目标。

长期建立的同质园活体收集对这一解释至关重要。采样树木已在同一位点生长约18–25年,减少了使地理分隔自然群体间比较复杂化的即时站点间环境对比。然而,该收集不等同于同步、随机化的种源试验:建植年份不同、来源信息不完整,且年龄、母体、发育和局部微环境效应无法正式分离。共享站点因此为采样种质间持久类群相关光谱分化提供证据,但不识别这些差异的因果基础或证明遗传分化。以往研究表明光谱表型可在作物和森林树木材料中恢复分子标记定义的遗传组,而同质园实验证明遗传身份、环境及其互作可在叶片反射率中联合检测。相比之下,本研究无分子数据或重复种源可用来划分这些变异来源。因此研究人员的推断限于广泛共享站点条件下持久的类群相关分子组成、散射相关和元素相关光谱差异。

模态特异性热图和树状图表明互补性在监督分类之前已显而易见。这对N. × leoni尤为清晰,即N. obliqua和N. glauca之间的天然杂交种:在NIR树状图中与N. glauca聚类(AU = 100,BP = 100),而在XRF树状图中与N. obliqua聚类(AU = 99,BP = 99)。该模式最好解释为跨模态与两个亲本物种的对比表型亲缘关系。这些亲缘关系是描述性的,可能源于生理、结构、元素稳态、种源相关或遗传过程。由于树状图是表型性的而非系统发育性的且此处未产生分子数据,该模式不能归因于差异亲本遗传或渐渗。更一般地,光谱相似性可能源于共享祖先、渐渗、叶片性状趋同、元素经济、叶习性或在具有 documented 杂交和网状进化历史的谱系中这些过程的组合。

超越这些终端亲缘关系,探索性分组分析表明两个光谱域在反映更广泛功能组织方面也存在差异。叶习性在XRF中比在NIR中结构更清晰,而亚属归属在两模态中均明显但XRF中再次更强。这种对比在生物学上合理,因为常绿和落叶策略在叶片寿命、养分投资和养分保留动态方面不同,这些性状构成更广泛叶片经济谱的一部分。衰老期间的养分再吸收同样与植物养分经济密切相关并在木本植物策略(包括常绿和落叶物种)间不同。由于XRF直接敏感于元素相关发射信号,养分获取、保留、再转运和分配的差异可能促成该光谱域中观察到的更强叶习性组织。近期杨属离子组学工作进一步证明叶片元素谱是具有实质性遗传组分的生物学结构化表型而非仅土壤化学的简单反映。

相反,NIR整合关于生化组成和散射相关属性的宽泛且部分重叠的信息,且性状-光谱关系在功能群和生物学背景间变化。因此NIR中较弱的叶习性分离并不意外,也不意味着叶习性在NIR表型中缺失;相反,其信号可能分布于多个协变化学属性而非通过更离散的元素相关谱表达。重要的是,叶习性和亚属对比均不能解释为独立效应,因为两属性在该收集中均由分类学身份结构化。因此其价值主要是解释性的:它们表明假山毛榉属类群间的光谱分化并非任意,而是部分与公认的功能策略和进化组织一致。

XRF结果在生物学上合理,因为离子组作为整合多元素系统运作而非一组独立浓度。矿物获取、转运蛋白活性、根-微生物互作、维管分配、叶片寿命和稳态调节共同塑造元素组成。在共享站点下,这些过程中的物种相关差异可产生持久多元素谱。Fe、Zn、Ca和Mn显示最大类群效应,而Mn、Fe、K和Si对缩减特征PLS-DA贡献最强。较弱且较不稳定的P响应也具有启发性:全类群集的统计显著性在平衡敏感性分析中未持续,强调将稳健元素分化与边缘发现分离的价值。

八信号模型性能仍远低于全XRF光谱。这一差异表明判别不限于主导整合峰。全光谱模型可利用峰形、肩部、次级线、连续谱结构和协调弱信号,这些在光谱压缩为少量峰面积时丢失。部分额外信息可能反映真实多元素模式,而部分可能源于与标准化粉末相关的基质和散射效应。因此全光谱应描述为强大的光谱指纹而非直接多元素浓度谱。峰面积经背景校正并在标准化制备下分析,但吸收、颗粒大小、堆积密度、厚度和几何形状仍可影响XRF响应。绝对浓度声明需要基质匹配校准或与基于ICP的测定比较。

NIR代表更连续和分布式的表型。约896–2515 nm的宽谱带整合许多协变化学成分和物理散射响应,这可降低近缘类群间的可分性,即使存在生物学意义变异。因此较低NIR准确率不应解释为分析失败。反射光谱已在栎属和种子植物中捕获遗传和系统发育结构,但信息信号强度和位置在谱系、环境和实验设计间变化。在本数据中,NIR以强自助支持解析若干终端关联并贡献足够非冗余信息以改变融合敏感性,但未匹配XRF实现的分类学判别。

O2PLS分解阐明为何两模态尽管预测准确性不同仍可互补。每区块近半数标准化变异由联合成分捕获(NIR 45.8%,XRF 52.1%),且配对联合得分显示中至强跨模态相关。由于物种标签未用于推导O2PLS潜在空间,共享表示中高类群相关变异(R2 = 0.605)提供独立描述性证据表明大部分NIR-XRF协变是生物学结构化的而非仅监督分类的结果。同时,NIR特异性和XRF特异性空间保留显著分类学组织(R2分别为0.212和0.300),表明互补性包含共同和模态特异性维度。低轮廓值警示不应将共享空间解释为九个离散簇;相反,类群沿重叠潜在梯度位移。同样,共享NIR和XRF载荷区域应解释为协变光谱域而非NIR直接检测单个元素的证据。

融合分析需要刻意谨慎的解释。特征缩放和区块均衡融合相对于XRF仅提高平衡分析准确率1.16和1.99个百分点,且经验分区水平范围包括零。区块均衡提高了N. glauca、N. macrocarpa和N. obliqua的平均敏感性但降低了N. dombeyi,而N. betuloides在所有模态中仍难以解析。因此融合模型改变了误差景观而非均一地改善每个类。区块均衡本身是分析加权选择:为标准化NIR和XRF区块分配相等预期总方差改变融合预测空间的几何形状并可重新分配类群间分类误差。N. dombeyi的降低敏感性说明了这一风险。此处略高的点估计和较低的重采样SD因此应视为数据集特异性而非区块均衡将普遍改善性能的证据。需要外部验证以确定小的内部优势是否在独立树木、收集、年份、仪器和站点间持续。

尽管平衡分析中每物种12株独立树木的重复对于捕获种内变异全貌可能显得适度,但研究设计旨在表征活体收集内连贯的类群水平分化而非估计全物种多样性。重要的是,生物学信号得到几条收敛证据线支持:模态特异性物种均值结构、强支持终端簇、重复交叉验证判别以及多个元素相关XRF信号间的实质性种间效应。这种收敛表明分类学分化不是单一分析程序的伪影。此类解释与新兴光谱生物学框架一致,后者将光谱视为能够连接分类学身份与跨生物尺度生化、结构和功能变异的整合表型。它也与近期同质园证据一致,表明叶片反射率在森林树木中保留可检测的遗传、环境和基因型-环境组分。因此,尽管种内重复适度,本研究提供生物学连贯且分析上可重复的共享站点条件下分类学和元素分化表征,同时避免关于全物种表型多样性的声明。

本研究中所有预测指标量化单一活体收集和一次采样活动内的内部判别;它们不是向其他园、年份、种源、仪器或野生群体可迁移性的估计。活体收集可提供卓越的分类学覆盖和长期研究基础设施,但其种质不一定捕获野外存在的种源结构或种内多样性。在本研究中,N. × leoni由六株独立树木代表,而三个N. nitida光谱源自一株生物学树木的分析重复;因此,平衡九物种分析提供预测性能的主要基准。重复种源的缺失和同步建植进一步妨碍类群、种源、年龄和持久发育效应的正式分离。采样也限于2025年1月8–9日,因此观察到的光谱模式代表一个南半球夏季物候期的叶片状态而不表征季节变异。最后,XRF特征代表标准化相对信号而非校准元素浓度。近期基于定量元素平台和数百基因型的杨属同质园离子组学研究已证明遗传对叶片元素组成的实质性贡献;研究人员的发现与生物学结构化元素表型相容但应视为比较微区XRF指纹而非定量离子组学。

综上所述,这些发现将配对NIR和微区XRF光谱分析定位为表征活体收集中互补类群相关变异的框架。在该收集内,全光谱XRF提供最强的独立类群判别,而NIR贡献互补分子组成和散射相关信息并在与XRF整合时产生温和额外增益。因此多模态采集的价值较少在于自动替代强单一模态而在于拓宽从同一标准化材料恢复的生物学维度。任何用于分类学验证、种质审计或保护决策的操作性使用将需要种源分辨采样和跨独立收集、季节和分析设置的外部验证。

研究结论部分翻译:配对NIR和微区XRF光谱分析在南美假山毛榉属长期迁地收集内恢复了类群相关变异的互补维度。全光谱XRF在收集内提供最强的独立判别,而NIR贡献非冗余信息且多模态整合仅在预测准确性上产生温和、类群依赖的增益。O2PLS进一步表明相当部分变异跨模态共享而模态特异性空间保留额外类群相关结构。这些发现支持光谱表型分析作为活体收集比较表征的工具框架,但它们未建立遗传分化或超出采样种质的一般分类学验证。在方法可用于种质审计或保护管理的操作性使用之前,需要独立收集、种源感知采样、重复季节测量和外部验证。
相关新闻
生物通微信公众号
微信
新浪微博
  • 搜索
  • 国际
  • 国内
  • 人物
  • 产业
  • 热点
  • 科普

热点排行

    今日动态 | 人才市场 | 新技术专栏 | 中国科学人 | 云展台 | BioHot | 云讲堂直播 | 会展中心 | 特价专栏 | 技术快讯 | 免费试用

    版权所有 生物通

    Copyright© eBiotrade.com, All Rights Reserved

    联系信箱:

    粤ICP备09063491号