整合转录组与DNA甲基化测序揭示圈养大熊猫雄性自然交配能力的分子调控机制及其保护意义

本研究通过整合血液转录组和全基因组DNA甲基化测序（WGBS）技术，系统揭示了圈养雄性大熊猫自然交配能力差异的分子基础。研究发现关键差异表达基因（DEGs）如ZPBP2、PIWIL4等参与GnRH信号通路、MAPK级联和嗅觉感知等过程，同时鉴定到启动子和基因体相关的差异甲基化区域（DMRs）与基因表达呈负相关。多组学整合分析表明，DNA甲基化通过调控轴突导向、半胱氨酸和蛋氨酸代谢、细胞自噬等繁殖相关通路影响交配行为。该研究为濒危物种繁殖行为的表观遗传调控提供了新见解，对优化圈养繁殖策略具有重要价值。

　　

1 引言

自然交配能力是圈养大型濒危哺乳动物维持繁殖能力和行为适应性的关键性状。与人工授精等辅助生殖技术相比，自然交配能更好地保留性选择、遗传多样性和行为完整性，并减少长期人工干预导致的行为异常和繁殖失败风险。在圈养环境中，自然交配行为的丧失被广泛认为是适应性下降的表现，恢复这种能力已成为迁地保护计划的核心焦点。尽管技术不断进步，许多大型濒危哺乳动物（包括大熊猫、北极熊、猎豹和云豹）在圈养环境中仍面临自然交配困难。例如，1998年至2018年间，圈养大熊猫的自然交配成功率普遍低于50%。约10%–30%的雄性大熊猫缺乏性动机或表现出较差的交配行为，这成为繁殖效率的主要障碍。虽然激素水平、精子质量、社会经验、行为表达、配偶选择、饮食和肠道微生物等多种外部因素与交配行为相关，但在受控条件下，雄性个体间交配能力的显著差异表明内在分子机制可能起决定性作用。

2 材料与方法

2.1 伦理声明

本研究中使用的血液样本是在年度常规健康检查期间采集的。所有大熊猫均未受到伤害，并在采样后恢复良好。研究方案获得北京师范大学动物福利与伦理委员会的批准（批准号：CLS-AWEC-E-2022-002-1）。

2.2 研究对象、血液样本采集和实验设计

成功与发情雌性完成爬跨、插入和射精，并使雌性产仔的雄性大熊猫被归类为有能力组。在交配季节多次尝试与若干发情雌性交配均失败，表现出对雌性不感兴趣、对雌性有攻击行为，或试图以错误姿势或方向爬跨的雄性被归类为无能力组。所有血液样本均在每年3月至5月的大熊猫交配季节采集。样本来自陕西省秦岭大熊猫研究中心的7只圈养大熊猫（4只有能力，3只无能力）、四川成都大熊猫繁育研究基地的8只大熊猫（3只有能力，5只无能力）以及中国大熊猫保护研究中心汶川卧龙神树坪基地的6只大熊猫（4只有能力，2只无能力）。根据样本可用性和年龄，最终选择18只个体（9只有能力，9只无能力），平均年龄为13.83±4.32岁，用于转录组测序分析。

2.3 RNA提取与文库构建

使用PAXgene Blood RNA Kit（Qiagen, Germany）提取RNA，并通过Agilent 2100 bioanalyzer（Agilent Technologies, CA, USA）评估RNA样品的纯度、浓度和完整性。文库构建以总RNA为起始材料，使用oligo(dT)磁珠富集带polyA尾的mRNA。随后将mRNA随机片段化，合成第一链cDNA，降解RNA链后合成第二链cDNA。纯化的双链cDNA经过末端修复、加A尾和连接测序接头，使用AMPure XP beads（Beckman Coulter, CA, USA）选择370–420 bp的cDNA片段，进行PCR扩增和纯化，最终获得文库。文库质量通过Qubit 2.0 Fluorometer（Thermo Fisher Scientific, MA, USA）定量和Agilent 2100 bioanalyzer检测插入片段大小进行评估。

2.4 转录组测序、质量控制、序列比对和表达分析

通过Illumina NovaSeq 6000平台进行测序，生成150 bp双端读长。原始测序数据经过过滤，去除带接头的读段、含模糊碱基（N）的读段和低质量读段（Phred质量值≤5的碱基超过50%）。清洁数据的Q20值达到98.10%±1.02%。使用HISAT2（v.2.0.5）软件将高质量序列比对到大熊猫参考基因组（NCBI accession ASM200744v3），平均93.27%的读段成功比对，其中88.89%为唯一比对。通过StringTie软件（v.1.3.3）进行转录本组装，使用featureCounts（v.1.5.0）计算每个基因的读段计数，并以FPKM（fragments per kilobase of exon model per million mapped fragments）表示基因表达水平。

2.5 基因功能分析和转录组差异表达分析

使用clusterProfiler软件（v.3.8.1）进行Gene Ontology（GO）功能富集分析和Kyoto Encyclopedia of Genes and Genomes（KEGG）通路富集分析。基于FPKM值计算样本间的Pearson相关系数（R2>0.8视为强相关）。主坐标分析（PCoA）和统计检验使用R软件（v4.1.0）的vegan包（v2.5–7）进行。差异表达基因（DEGs）的筛选标准为|log₂(FoldChange)|>1且错误发现率（FDR）调整后p值≤0.05。随机森林分析使用R包randomForest（v.4.7–1.2）进行，以识别区分有能力组和无能力组的关键特征基因。

2.6 甲基化样本制备、测序数据获取和预处理

每个样本取100 ng基因组DNA，与0.5 ng未甲基化lambda DNA混合，使用Covaris S220超声破碎仪片段化至200–400 bp。使用EZ DNA Methylation-Gold Kit进行未甲基化胞嘧啶的亚硫酸盐转化。文库构建后，通过Agilent 5400系统和qPCR评估质量，合格的文库在Illumina NovaSeq 6000平台上进行双端测序。原始测序数据经过FastQC（v0.11.8）评估和质量过滤（使用fastp v0.23.1），生成高质量清洁数据。使用Bismark（v0.24.0）将清洁的亚硫酸盐处理读段比对到参考基因组，并去除PCR重复。通过bismark_methylation_extractor（--no_overlap）提取胞嘧啶位点的甲基化信息，并根据lambda DNA中未甲基化胞嘧啶估计亚硫酸盐转化效率。

2.7 甲基化水平和差异甲基化分析

使用二项检验评估每个胞嘧啶位点的甲基化显著性，错误发现率（FDR）调整后p值<0.05的位点定义为甲基化。为评估全基因组甲基化水平，将基因组划分为10 kb窗口，计算每个窗口内甲基化位点的甲基化水平（ML）。使用DSS软件（v2.12.0）识别差异甲基化区域（DMRs），DMRs根据其基因组位置分为基因编码区（转录起始位点TSS到转录终止位点TES）和启动子区（TSS上游2 kb）。使用GOseq进行GO富集分析，使用KOBAS软件（v3.0）进行KEGG通路富集分析。使用Homer（findMotifsGenome.pl v4.11）进行motif分析，以识别基因体和启动子区中长度8、10、12和14核苷酸的前25个显著motif。

2.8 DNA甲基化与转录组的整合分析

在整合全基因组DNA甲基化和转录组数据的基础上，采用多种可视化和统计方法系统评估它们之间的关联。首先，在染色体水平上，合并每组内的生物学重复样本，计算CG、CHG和CHH上下文的甲基化密度以及转录组测序读段密度，使用Circos图（v0.62-1）展示跨染色体的空间分布。其次，在基因及侧翼区域水平，根据FPKM值将所有基因分为四类：无表达（<1）、低表达（1至第一四分位数）、中表达（第一至第三四分位数）和高表达（>第三四分位数）。在基因体及其±2 kb侧翼区域内，各分为50个等分区间，计算每个区间内胞嘧啶的平均甲基化水平，并绘制不同表达水平下的甲基化模式图。第三，根据启动子（TSS上游2 kb）和基因体中的甲基化水平将基因分为五个分位数（20%、40%、60%、80%和100%），比较组间的表达分布（log₂(FPKM+1)），并使用箱线图和散点图分析甲基化与基因表达的相关性。对于DMR相关基因分析，通过Wald检验（阈值0.2，p.threshold=1e-5）识别DMRs，并将其相关基因集与DEGs（|log₂FC|≥1）取交集。对重叠基因进行层次聚类热图和线图分析，展示跨基因体及±2 kb侧翼区域的甲基化动态。

2.9 统计分析

为识别有能力组和无能力组之间血液转录组中的显著DEGs，首先计算两组间基因表达的log₂转换折叠变化（log₂FoldChange），筛选|log₂FoldChange|>1的基因进行表达差异分析。使用Benjamini–Hochberg（BH）方法调整p值以控制错误发现率（FDR）。基于DEG分析结果生成火山图。随机森林分析使用R包randomForest（v.4.7–1.2）基于基因FPKM值进行。对DEGs进行Pearson相关分析以探索共表达模式，使用ggplot2包（v.3.5.0）可视化结果。基于Bray–Curtis相异度指标使用PERMANOVA评估组间转录组数据差异，通过vegan包（v2.6–8）的adonis2函数实现，显著性基于999次置换检验。使用DSS软件（DSS-single模型）识别DMRs，该模型通过平滑相邻CpG位点之间的空间相关性、整合测序读长深度以及使用β-二项分布建模组内变异，来提高单点甲基化估计和DMR检测的准确性和敏感性。

3 结果

3.1 与雄性大熊猫自然交配能力相关的基因表达和差异表达基因（DEGs）

在圈养雄性大熊猫的血液转录组中共鉴定到33,286个基因，包括1954个新基因，用于基因表达量化。有能力组和无能力组样本的FPKM值相关性分析显示组内相关性强（R2>0.9）。基于Bray–Curtis相异度的主坐标分析（PCoA）显示两组间无显著分离（F=1.52，R2=0.087，p=0.166；PERMANOVA）。当分别分析四川和秦岭两个繁殖基地的样本时，得到一致结果。尽管如此，使用|log₂(FoldChange)|>1且调整后p值≤0.05的阈值，在有能力组和无能力组之间鉴定到91个DEGs。其中，32个基因在有能力组中显著下调，59个上调，包括关键基因如LOC117800792、ZPBP2、DNAH7和PIWIL4。随机森林分析基于平均减少精度值识别出前10个对组分类贡献最大的特征，不包括假基因，显著贡献者包括LOC117795464、LOC100483191、ICA1L、LOC105234470、LOC117795679、LOC117797362和LOC117795847。基于DEGs的共表达分析显示五个枢纽基因具有最高度值：LOC105235995、RAB3C、LOC117800801、ZPBP2和SOAT2。有能力组中上调基因的GO富集分析突出了一些与生殖和信号功能相关的术语，包括生殖过程、繁殖、多生物过程、信号转导活性和信号传导，表明与生殖和信号过程的关联。更广泛的生物过程如对刺激的反应、代谢过程和细胞过程也被富集，表明生理反应性和代谢活性升高。值得注意的是，ZPBP2（透明带结合蛋白2）是一个在生殖GO术语中注释的基因，编码一种参与结合卵子透明带的精子蛋白，并包含一个保守的PF07354结构域，对精子-卵子识别至关重要。DEGs的KEGG通路分析进一步揭示了与繁殖和类固醇生成相关通路的富集，包括类固醇激素生物合成、卵巢类固醇生成和胆固醇代谢。Wnt信号通路和Notch信号通路的富集也表明可能涉及睾丸发育和精子发生。通过比对NR数据库，我们共鉴定到320个与雄性大熊猫交配行为直接相关的基因，涵盖六个功能类别：雄激素信号、GnRH信号、精子质量（包括精子发生、运动、结合和受精）、嗅觉受体、多巴胺信号和生殖细胞发育。其中，GnRH通路相关基因的表达水平（FPKM）排名前20。蛋白质-蛋白质相互作用网络分析进一步突出了10个核心蛋白——RAF1、MAP2K1/2、MAP3K1、RAS家族、SOS1和PIK3CA/D，表明它们在调节交配行为的关键信号通路中起核心作用。

3.2 WGBS分析

对四川成都大熊猫繁育研究基地的雄性大熊猫进行了全基因组亚硫酸盐测序（WGBS）。测序数据质量高，原始读段数约236–275百万，过滤后清洁读段比例在86.85%到87.92%之间。质量分数高（Q20>96.7%，Q30>90.3%），GC含量稳定（21.3%–22.1%），亚硫酸盐转化率始终高于99.6%。所有样本的CG甲基化水平和密度高且稳定，而CHG和CHH甲基化水平和密度显著较低，并且在样本间和基因组窗口间变异性更大。在三种序列上下文中，甲基化主要集中在CG位点，CG甲基化密度高达99.17%，而CHG和CHH上下文的最大密度分别仅为0.25%和0.21%，反映了哺乳动物基因组的典型甲基化模式。在基因组功能区域内，外显子和内含子等基因体区域的CG甲基化水平显著较高，而启动子和非翻译区（UTR）区域表现出显著的低甲基化。尽管CHG和CHH甲基化水平全局较低，但在启动子和5' UTR区域观察到两组间存在轻微差异。CG甲基化表现出典型的双峰分布模式，即转录起始位点（TSS）附近甲基化降低，基因体内水平升高，随后在转录终止位点（TES）附近下降。值得注意的是，有能力组基因区域上游和下游的CHG和CHH甲基化水平始终高于无能力组。

3.3 DMR分析

在DMRs分析中，所有样本的CG甲基化水平表现出高度一致性（Pearson's R2>0.8），而CHG和CHH上下文中的甲基化水平相关性低，R2值分别不超过0.1和0.12。有能力组和无能力组之间的甲基化差异主要出现在CG序列上下文，其占总DNA甲基化的80%以上。两组之间共鉴定到6656个DMRs，包括6321个CG型DMRs。然而，CG、CHG和CHH上下文的整体甲基化水平分布在组间未显示显著差异。在CG上下文中，大多数DMR相关基因锚定在基因体区域（n=2930），而只有五个基因同时在所有三种上下文（CG、CHG和CHH）中被注释。类似地，在启动子区域，CG上下文也产生了最多数量的DMR相关基因（n=979），其中31个基因在所有三种序列上下文中重叠。基于CG上下文DMR甲基化水平的聚类热图清晰地区分了有能力和无能力个体，而CHG和CHH上下文由于其全局低甲基化水平，未揭示明显的组间差异。在所有鉴定到的DMRs中，76个区域的组间绝对甲基化差异大于0.5，包括28个高甲基化区域（全部在CG上下文）和48个低甲基化区域（主要在CG，2个在CHG）存在于有能力组。总共有3045个基因在各种基因元件（如外显子、内含子和UTR）中映射到DMRs，990个基因与位于启动子区域的DMRs相关。

3.4 DMR相关基因的功能富集分析

基于DMRs与基因组功能区域的重叠，我们对基因体区域DMRs相关基因进行了KEGG通路富集分析。在CG序列上下文中，DMR相关基因显著富集（p<0.05）于几个信号通路，包括轴突导向、钙信号通路、谷氨酸能突触、Rap1信号通路、Ras信号通路、Wnt信号通路和肌动蛋白细胞骨架的调节。进一步分层分析显示，在有能力组中，高甲基化（Hyper）DMR基因主要富集于轴突导向、肌动蛋白细胞骨架的调节和Rap1信号通路，而低甲基化（Hypo）DMR基因显著富集于癌症通路、多巴胺能突触、胆碱能突触、昼夜节律授时、Ras信号通路、磷脂酶D信号通路和钙信号通路。在CHH序列上下文中，DMR相关基因显著富集于氧化磷酸化、帕金森病、产热、逆行内源性大麻素信号、代谢通路、阿尔茨海默病、亨廷顿病、非酒精性脂肪肝病（NAFLD）和心肌收缩。对于CHG上下文，富集通路包括逆行内源性大麻素信号、氧化磷酸化、帕金森病和产热。此外，具有启动子区域DMRs的基因显著富集于MAPK信号通路、肌动蛋白细胞骨架的调节、Rap1信号通路和cAMP信号通路。

3.5 DNA甲基化与转录组谱的整合分析

总体而言，DNA甲基化水平和基因密度的染色体分布模式在有能力组和无能力组之间具有可比性。然而，在X染色体（NC_048238.1）上观察到CG甲基化密度存在显著差异。通过整合甲基化组和转录组数据集，在DEGs和DMR相关基因之间共鉴定到3094个重叠基因。其中，47个基因表现出甲基化增加和表达减少的负相关关系，而45个基因表现出甲基化减少和表达增加的关系，基于0.2的差异甲基化阈值和|log₂(Fold Change)|>1的阈值。在基因体区域，70个和96个基因分别表现出DMRs与DEGs之间的正相关和负相关，而在启动子区域，21个和34个基因显示出这种关联。CG上下文DMRs与基因表达的相关性和聚类分析显示，组间启动子甲基化差异比基因体区域观察到的差异更为明显。具有基因体甲基化的重叠基因的功能富集表明显著的GO术语，如“蛋白质代谢过程”、“细胞生物合成过程”和“含氮化合物代谢过程的调节”，以及KEGG通路包括轴突导向、趋化因子信号通路、鞘脂信号通路、凋亡/自噬、精氨酸和脯氨酸代谢以及O-聚糖生物合成。对于具有启动子相关DMRs的基因，GO富集见于“蛋白质代谢过程”、“核苷酸磷酸结合”和“生物合成过程的调节”，而KEGG通路富集于糖胺聚糖生物合成、C型凝集素受体信号、溶酶体、轴突导向、Hedgehog信号、半胱氨酸和蛋氨酸代谢、肌醇磷酸信号和VEGF信号。

4 讨论

本研究整合转录组和WGBS数据，研究了圈养雄性大熊猫自然交配能力的分子基础。具体而言，我们的研究发现，尽管全局基因表达谱在有能力组和无能力组之间没有显示出显著分离，但差异表达分析揭示了一个与精子发生、精子-卵子识别和睾丸发育相关的基因子集（例如ZPBP2、PIWIL4和DNAH7）。这表明自然交配能力可能受少数关键功能基因的调控。此外，WGBS分析展示了典型的哺乳动物DNA甲基化模式，即高且稳定的CG甲基化和低的非CG甲基化（CHG/CHH），但揭示了组水平的细微差异，特别是在基因侧翼区域的CHG和CHH上下文中。重要的是，DMRs主要出现在CG上下文中，并且显著富集于参与轴突导向、生殖信号和激素相关通路的基因。甲基化和基因表达的整合分析进一步鉴定出表现出经典甲基化-表达反相关模式的基因，特别是在启动子区域内，突出了转录活性的潜在表观遗传调控。此外，在X染色体上观察到的差异甲基化暗示了性染色体连锁调控元件的可能参与。总之，这些发现表明，雄性大熊猫的自然交配能力可能受到与生殖功能、神经信号和细胞通讯相关的基因的协同转录和表观遗传调控的影响。

我们的圈养雄性大熊猫具有不同自然交配能力的转录组分析揭示了一个与生殖生理、神经内分泌信号和细胞调控密切相关的DEGs子集。值得注意的是，编码透明带结合蛋白2的基因ZPBP2在有能力组中显著上调。这种蛋白质对精子-卵子识别至关重要，并且已被报道为多种哺乳动物精子顶体完整性和生育力的标志物。它的上调表明具有成功交配行为的雄性精子功能能力增强。此外，共表达和基于机器学习