梅花鹿粪便中线粒体DNAD-loop区异质性的研究

A Study on the Heteroplasmy in D-loop Region of Sika Deer Fecal Mitochondrial DNA

摘要

线粒体DNA（mtDNA）的异质性（heteroplasmy）是多种不同序列的mtDNA拷贝共存于一个个体的现象，对线粒体、细胞和组织的机能产生明显的影响。粪便是野生动物研究中最易得的材料，含有动物肠上皮脱落的细胞及其DNA，其中的mtDNA异质性可以反映个体在各种因素影响下肠上皮细胞DNA的稳定性，为揭示个体的适合度、生态对策等提供新的线索。梅花鹿（Cervus nippon）是森林生态系统中的关键物种，其肠上皮中mtDNA异质性的特点尚未见报道。本研究用27只6～96月龄的圈养梅花鹿粪便样本，采用PCR扩增子高通量测序技术，对mtDNA D-loop区的异质性进行定量分析。结果显示：D-loop区广泛分布着异质性位点，但相对集中于第175～368位碱基的高变区（HyperHPR），而第597～687位碱基区间则较少发生，为保守区（HypoHPR）。HyperHPR异质性突变的多样性（HDα）极显著高于HypoHPR（P < 0.001）。转换（T_S）是主要碱基替换类型，其发生频率远高于颠换（T_V）。HyperHPR中T_S/T_V平均为0.11，而HypoHPR平均为0.19，说明异质性突变在HypoHPR受到更强的选择性清除。无论是在D-loop区、HyperHPR与HypoHPR，雌雄个体之间异质性突变指标均无显著差异（D-loop区，\left|t\right| = 0.273，P = 0.787；HyperHPR，\left|t\right| = 0.409，P = 0.686；HypoHPR，\left|t\right| = 0.793，P = 0.435）。此外，HyperHPR的异质性频率与月龄之间存在显著正相关关系（R² = 0.132，P = 0.035），使整个D-loop区也具有近乎显著的趋势（R² = 0.104，P = 0.056）。但是，异质性突变各项指标均与性别无关（P = 0.435～0.787）。以上结果表明，肠上皮细胞的D-loop区作为接近中性进化的非编码区，能够保留大量的异质性突变，但对T_V突变产生明显的选择性清除，特别是保守区域对T_S和T_V均有强烈的清除作用，从而使异质性突变形成的mtDNA多样性控制在一个相对正常的范围内。年龄是影响异质性水平的重要因素，性别的影响较小。研究结果为利用粪便研究肠上皮的稳定性与环境的关系提供了新的视角。

Abstract

The heteroplasmy of mitochondrial DNA （mtDNA） is a phenomenon that multiple variants coexist in an individual. It affects the functions of mitochondria， cells， and tissues. Feces is the most conveniently collectable material from wildlife， containing exfoliated gut epithelial cell and their DNAs. It could be expected that the heteroplasmy of fecal mtDNA could be an indicator of genetic stability of epithelium under various environmental factors， providing new clues into individual fitness and effectiveness of ecological strategies. The sika deer （Cervus nippon） is a key species in forest ecosystems. The gut epithelial heteroplasmy has not yet been studied. The present study characterized the heteroplasmy in D-loop region of fecal mtDNA using fecal samples of 27 captive sika deer between 6 and 96 months old， and high-through put sequencing of PCR products. The results showed heteroplasmic mutation wide spread throughout the D-loop， but concentrated in the hyperheteroplasmic region （HyperHPR） ranging from the 175th to 368th nucleotide. The region between the 597th and 687th nucleotide was a hypoheteroplasmic region （HypoHPR） where carried very few mutations. The α diversity of heteroplasmy diversity （HDα） of HyperHPR was extremely significantly higher than that of HypoHPR （P < 0.001）. Transition （T_S） was the preponderant type of nucleotide substitution， of which the occurrence was significantly higher than transversion （T_V）. The ratio T_S/T_Vwas averaged 0.11 in HypoHPR and 0.19 in HyperHPR， suggesting the mutations were more selectively removed from HypoHPR. There was no significant difference between males and females for the entire D-loop region （\left|\mathit{t}\right| = 0.273， P = 0.787） and the HyperHPR （\left|\mathit{t}\right| = 0.409， P = 0.686） and HypoHPR （\left|\mathit{t}\right| = 0.793， P = 0.435）. Meanwhile， the frequency of heteroplasmy in HyperHPR was significantly correlated with age in month （R² = 0.132， P = 0.035）， leading to a marginal significance in D-loop region （R² = 0.104， P = 0.056）. However， no relatedness was observed between sex and heteroplasmic indices （P = 0.435 - 0.787）. The above results suggest that D-loop region of gut epithelium could accommodate a large quantity of heteroplasmic mutations though it generally evolves neutrally. However， selective elimination acts strongly for T_V across the whole region， and even stronger for both T_S and T_V in HypoHPR， which results in the diversity of heteroplamy laying below an acceptable threshold. Age is an influencing factor of heteroplasmic mutation， while sex is not. In general， this study provides a novel angle to regard the relationship between gut epithelium stability and environmental factors.

Graphical abstract

关键词

梅花鹿 / 粪便DNA / 线粒体异质性 / 年龄

Key words

Sika deer（Cervus nippon） / Fecal DNA / MtDNA heteroplasmy / Age

线粒体是真核细胞能量代谢的枢纽^［1］，负责能量转换和ATP的生产，同时还参与调节细胞内环境^［2］、促进信号传递和控制细胞的生长与凋亡^［3］，以及参与免疫应答过程^［4］。线粒体拥有独立环状基因组，称为线粒体DNA（mitochondrial DNA，mtDNA）^［5］，位于线粒体基质中。在哺乳动物中，mtDNA长约16.5 kb，包含37个基因，编码13个电子传递链的蛋白亚基、22种转运RNA（tRNAs）和2种核糖体RNA（rRNA）^［6］。

一个哺乳动物的细胞内含有数百到上千个mtDNA的拷贝，这种多拷贝属性赋予了细胞对能量需求的快速响应能力^［7］。与核基因组相比，mtDNA表现出更高的突变率和有限的修复机制。由于线粒体始终处于更新中，mtDNA的突变可以随着其增殖和消除而被扩增或消除，因此，随着时间的推移，细胞、组织和器官中出现序列不同于祖代mtDNA的拷贝，这种现象被称为mtDNA异质性（heteroplasmy）^［8］。每个异质性突变或多或少会影响线粒体以及细胞和组织的功能，进而影响物种对环境变化的适应能力，甚至引发严重的疾病^［9-10］。2023年，Wang et al.^［11］把细胞、组织、器官、系统乃至个体所包含的mtDNA突变的总和看成一个总体，提出了线粒体基因库（mitogene pool）的概念，指出其大小（即多样性的总量）可以反映各个水平上突变和突变清除的活跃程度，预示着对线粒体、细胞和个体整体机能的变化潜能，在生理学、生态学、进化生物学以及保护生物学中具有巨大的应用潜力，具有进化意义。描述线粒体基因库的特点，度量其大小，研究其与动物体某些性状的关系，可以从一个全新的视角来重新看待发育、衰老、自然选择、个体适合度和性状表型等经典生物学问题。

老年动物肠道的结构和功能往往呈现出退行性变化，如肠壁和黏膜层变薄、分泌机能下降等^［12-13］。其mtDNA的异质性随着年龄的增加具有增加的趋势，主要表现为特定区域突变位点数、突变的碱基组成和有害突变等的变化。虽然众多mtDNA拷贝中突变的碱基占比很小，但是从其组成来说，通常显示出强烈的负选择信号，即只保留影响效应较小的转换（transition，T_S），而清除了大量的颠换（transversion，T_V）^［14-15］。这些趋势预示着发育和衰老等进程中所伴随的细胞能量转换机能的微妙变化。但是，老年动物肠道的结构和功能这一变化趋势并不是简单的组织耗损，而是在细胞、结缔组织等结构性更新的进程中逐渐产生的替换能力不及损耗速度所致。那么，线粒体基因库如果不能维持一个健康的状态，也会直接影响肠道机能，使老年动物更脆弱，反过来会加剧衰退的进程。相反，如果线粒体基因库可以维持一个较好的状态，则会缓冲结构性退化带来的负面影响。

粪便是野生动物研究中最常见的非损伤性材料，特别是对野外难以接近的动物，几乎是唯一的生物学材料。粪便中含有较多的肠上皮脱落细胞^［16］。这些细胞与其他体细胞一样，含有大量的mtDNA拷贝，同样也可能发生异质性突变。解析粪便mtDNA的异质性，就可以揭示肠道上皮组织线粒体基因库的特征，进而与年龄、营养和毒素等因子进行关联分析。虽然粪便中的宿主DNA经过富集后可以用于开展高通量测序^［17］，但是粪便DNA的异质性解析仍然报道较少，缺少成熟的经验。

野生梅花鹿（Cervus nippon）为国家一级重点保护野生动物^［18］。在野外环境中，相比其他材料样本，粪便更易得，研究粪便中mtDNA的异质性对揭示其生物学特性具有重要意义。本研究以不同年龄梅花鹿的粪便DNA为材料，通过测定变异位点的数量、密度和频率等传统的异质性指标^［19］和新近提出的基于α多样性的新指标^［20］，系统量化了D-loop区在个体水平的基因库特点，初步探讨这些特点与年龄的相关性。

1 材料与方法

1.1 样本采集

从梅花鹿养殖场中采集27只6～96月龄的梅花鹿粪便样本，包括雌性15只，雄性12只。所有样品均为新鲜粪便，用一次性PE手套直接抓取粪团，放入自封口塑料样品袋中，立即转入冰箱中冷冻保存，统一由冷链运输至实验室，再转到冰箱中保存备用。操作时一个样本更换一次手套，防止交叉污染。

1.2 粪便DNA提取

使用无菌手术刀从粪便表层切取150～200 mg的样品，置于2 mL无菌离心管中。采用本团队创建的粪便中富集宿主DNA的PEERS法进行总DNA的提取。主要步骤：向装有粪便样本的离心管中加入500 μL 10 mmol/L的PBS溶液，充分浸泡约10 min，加入125 μL 5%的SDS，震荡，使粪便与溶液快速混合均匀，瞬时离心后于室温下静置5 min，3 000 r/min离心10 min，弃底部沉淀，将全部上清液转移至新的1.5 mL离心管内。用AxyPrep Multisource Genomic DNA Miniprep Kit试剂盒（Axygen，美国）从上清中提取总DNA。用NanoDrop 2000超微量紫外分光光度计（Thermo Fisher，美国）测定总DNA浓度。

1.3 D-loop区mtDNA异质性分析

1.3.1 D-loop区段的PCR扩增

以NCBI_GenBank中梅花鹿mtDNA序列为参考（OP_834097.1）设计4对PCR引物（表1），分别扩增其D-loop区内的4个不同片段，经拼接后可覆盖完整的D-loop区。引物序列经过NCBI-BLAST排除扩增核拷贝（Numts）的可能。各个片段的扩增均在50 μL体系中进行，包含2 × Rapid Taq Master Mix（南京诺唯赞医疗科技有限公司，中国）25 μL、上游和下游引物（10 μmol/L）各2 μL、DNA模板5 μL（~ 100 ng总DNA）和去离子水16 μL。用Eppendorf Nexus GSX1扩增仪（Eppendorf，德国）进行扩增，反应程序：95 ℃预变性5 min；94 ℃变性30 s，表1中相应温度下退火30 s，72 ℃延伸30 s，30个循环；最后72 ℃延伸10 min。扩增结束后，取5 μL扩增产物在2.0%琼脂糖凝胶上，100 V电压下电泳15 min，在紫外光下观察扩增产物的分布情况。

1.3.2 PCR产物的高通量测序

对每份样品的扩增产物进行切胶回收并纯化（赛文创新（北京）生物科技有限公司），用MGIEasy酶切DNA文库制备试剂套装（华大智造，中国）对待测DNA样品进行WGS测序文库构建。简要步骤为：取ERAT Buffer和ERAT Enzyme Mix加入到反应液中对DNA进行末端修复和添加dA尾；在反应产物中加入文库制备试剂盒中的Ligation Buffer与DNA Ligase，使之与DNA Adapters连接，然后使用DNA Clean Beads将连接产物进行磁珠纯化回收；用PCR Enzyme Mix加上PCR Primer Mix配置而成的PCR反应液，在50 μL体系内，按照前述反应程序扩增3个循环。全部流程均按照试剂盒说明书操作。

完成以上步骤后，通过毛细管电泳（Aligent 5200，美国）检测建库产物中主要目的片段的分布情况，使用双链荧光定量法，按照Qubit^® dsDNA HS Assay Kit荧光定量试剂盒的操作说明书对PCR纯化后产物进行定量，要求混合后DNA样本总量不少于 1 pmol，总体积不超过48 μL。用MGIEasy环化试剂盒（华大智造，中国）对DNA文库进行处理，该过程先后经历变性、单链环化和酶切消化3个步骤。将环化产物进行磁珠纯化回收，最后按照Qubit^® ssDNA Assay Kit的操作说明对酶切消化后产物定量，酶切消化产物产量（ssDNA） / PCR投入量（dsDNA）不小于7%。

最后，通过MGIDL-200H便携加样器（华大智造，中国）加载DNB，再由MGISEQ 2000 RS测序平台（华大智造，中国）进行PE150 bp（双端）测序。使用Fastp工具对原始下机数据（FASTQ格式文件）进行初步质量评估，去除含有接头、低质量核苷酸和不可识别核苷酸（N）的reads，保存clean reads备用。

1.3.3 mtDNA D-loop区异质性特征分析

用Trimmomatic软件通过命令行脚本对测序产生的双端序列数据进行质控处理，滤除低质量序列，切除序列接头。使用生物信息学软件MEGAHIT对每个样品的序列数据进行组装，生成一个高质量的D-loop区基因组序列。然后，利用BWA的mem算法将所有测序reads片段比对到组装好的梅花鹿线粒体内参基因组上。在比对过程中，通过设置兼容GATK的参数优化分析效率，输出文件由SAM格式转换为更高效的BAM格式。用Samtools对BAM文件进行排序、重复标记和重复去除，以减少假阳性变异的检出。用GATK的Mutect2工具在“mitochondria-mode”模式下分析线粒体DNA（chrM）的变异。通过指定输入输出文件和相关参数进行变异信息分析，并用FilterMutectCalls工具过滤假阳性变异，最终保留“FILTER”列为“PASS”的变异，具体参数为：（1）最大错误发现率设置为0.05；（2）“Alt”列读取的碱基质量设置为20，映射质量设置为30。

1.3.4 统计分析

在进行异质性分析时，根据Nei et al.^［21］提出的模型对D-loop区整体进行划分，从中定义异质性高变区（hyperheteroplasmic region，HyperHPR）与异质性保守区（hypoheteroplasmic region，HypoHPR）。HyperHPR的定义为D-loop区内多个连续异质性位点的平均频率高于D-loop区平均值的区域；HypoHPR采用更严格的标准，定义为多个连续异质性位点的平均频率显著低于0.05的区域。按照D-loop区整体以及HyperHPR、HypoHPR探究异质性的分布特点，包括异质性位点数量、异质性位点的密度、各位点的异质性频率和各突变类型频率，各指标的定义及计算方法如下。

（1）异质性位点的数量，异质性位点是指在同一DNA分子的特定位置（或同一细胞的不同DNA分子中，mtDNA），存在两种或多种不同的等位基因的现象。统计存在不同类型核苷酸的多态性位点数量。

（2）异质性位点的密度，指每100个碱基对内的异质性位点的数量。

（3）异质性频率，在特定位点上每种异质性碱基突变的发生频率，用每个异质性碱基的read数在该位点上总read数中的占比来衡量。

（4）碱基突变类型，包含T_S、T_V、片段插入（Insertion，I_ns）和缺失（Deletion，D_el）的数量。

（5）T_S/T_V，同一个碱基位点上发生的T_V突变数与T_S突变数的比值。

此外，整理27只梅花鹿的测序数据，采用Liu et al.^［20］的方法，希尔数q = 2，计算每只个体D-loop区整体、HyperHPR和HypoHPR异质性的α多样性指数（heteroplasmy diversity α，HDα）。通过Welch’s方差检验3个区域的α多样性指数是否存在显著差异。最后，将上述各项指数与月龄进行线性回归，检验二者之间的相关性，显著水平α = 0.05。

2 结果

2.1 mtDNA D-loop区异质性空间分布格局

从27只梅花鹿粪便中提取的总DNA浓度均在10～20 ng/μL，4对引物（Mt-D1、Mt-D2、Mt-D3和Mt-D4）的扩增产物在2.0%琼脂糖凝胶上均显示为单一而明亮的条带，大小符合设计预期。以每份样品的PCR产物经高通量测序后组装的D-loop参考序列共995 bp，与GenBank中梅花鹿的D-loop序列完全一致。将每份样品的测序数据分别与其比对，从27只梅花鹿中共得到异质性位点108个。从频率分布来看，异质性位点广泛分布在整个D-loop区内（图1）。D-loop区上确定1个HyperHPR区域和1个HypoHPR区域，HyperHPR位于第175～368位碱基区间；HypoHPR位于第597～687位碱基区间。相对于HyperHPR，HypoHPR的异质性位点密度较低，异质性位点的数量与每位点异质性频率的平均值明显低于前者（表2）。特别值得注意的是，D-loop区内普遍存在一个特殊的高频率插入/缺失位点，即第837位碱基，其类型为T/TC/TCC/TCCC/TCCCC。

2.2 mtDNA D-loop区异质性碱基替换特点

从D-loop区整体来看，所有个体的异质性突变中T_S占绝对优势，而T_V的频率极低（表3）。具体来说，共鉴定出203次碱基替换事件，其中197次为T_S，仅6次为T_V。个体上检测到的T_S为1～20个，平均为（7.30 ± 5.34）个；T_V为0～2个，平均为（0.22 ± 0.51）个；T_S/T_V为2～18，平均为（1.61 ± 4.41）。超过80%的个体中只发生T_S事件而无T_V事件。此外，D-loop区还观察到大量的I_ns、D_el事件，共计观察到I_ns事件106次，平均每个体（3.93 ± 2.57）次；D_el事件52次，平均每个体（1.92 ± 0.87）次。这些事件几乎都发生在HyperHPR和HypoHPR之外的区域，HyperHPR中I_ns与D_el事件发生的概率很低，仅为（0.07 ± 0.27），而HypoHPR甚至不存在片段的I_ns与D_el事件。另一方面，HyperHPR中T_S/T_V为0～3，平均值为0.11；而HypoHPR为0～5，平均值为0.19，说明HypoHPR受到更强的选择压力。

2.3 D-loop区异质性多样性

结果如图2，不同个体D-loop区总体的HDα为3.8～23.4，平均为（10.5 ± 6.7）；HyperHPR的HDα为0～19.0，平均为（2.4 ± 4.4）；HypoHPR的HDα为0～5.9，平均为（0.4 ± 1.2）。Welch’s方差检验结果表明，不同区域在HDα水平存在显著差异（P < 0.001）。

按照性别分组，统计出雄性个体D-loop区总体的HDα为4.0～23.0，平均为（4.8 ± 6.6），雌性个体D-loop区总体的HDα为3.8～23.4，平均为（5.7 ± 6.9）；其中雄性个体HyperHPR的HDα为0～19.0，平均为（1.3 ± 5.3），雌性个体HyperHPR的HDα为0～13.0，平均为（1.2 ± 3.7）；雄性个体HypoHPR的HDα为0～1.0，平均为（0.1 ± 0.4），雌性个体HypoHPR的HDα为0～5.9，平均为（0.3 ± 1.5）。无论是在D-loop区、HyperHPR与HypoHPR，雌雄个体之间均无显著差异（D-loop区，\left|t\right| = 0.273，P = 0.787；HyperHPR，\left|t\right| = 0.409，P = 0.686；HypoHPR，\left|t\right| = 0.793，P = 0.435）。

2.4 mtDNA D-loop区异质性与年龄的关系

把每只个体的异质性位点的密度、异质性频率、T_S、T_V、I_ns和D_el数、T_S/T_V以及HDα与相应的月龄进行线性回归分析，发现HyperHPR的异质性频率与年龄呈显著的正相关（R²= 0.132，P = 0.035），D-loop的异质性频率与年龄接近显著相关（R² = 0.104，P = 0.056，图3）。此外，整个D-loop区、HyperHPR以及HypoHPR的其他所有指数与月龄之间均没有显著的相关性（D-loop区：R² = 0～0.064，P = 0.107～0.615；HyperHPR：R² = 0～0.092，P = 0.068～0.788；HypoHPR：R² = 0～0.012，P = 0.261 ~0.606）。

3 讨论

D-loop区是mtDNA上几乎接近中性进化的区域，比其他区域有着更为丰富的突变^［22］，这使其成为探索线粒体异质性的理想区域。因为有着广泛的异质性，由扩增子高通量测序数据组装的D-loop区的序列是每个个体众多单倍型的一致性序列（consensus），以每个个体自己的一致性序列作为参考，在检测异质性突变时可在最大程度上减少假阳性信号。

从景观角度，D-loop区广泛分布着异质性位点，尤其是HyperHPR（第175～368位碱基），展现了异质性频率的高值。但是，D-loop区上也存在明显的异质性的低谷区域HypoHPR，较为典型的是第597～687位碱基，突变及其频率都较低（图1）。从碱基替换类型来看，D-loop区的T_S事件发生频率远高于T_V事件，尤其是HypoHPR，仅从当前的样本中观察到1次T_V事件。此外，D-loop区还观察到I_ns事件共计106次，平均每个体（3.93 ± 2.57）次；D_el事件52次，平均每个体（1.92 ± 0.87）次。由表3可知，HyperHPR中I_ns与D_el事件发生的概率很低，仅有（0.07 ± 0.27），而HypoHPR甚至不存在I_ns与D_el事件。HyperHPR区域内突变位点碱基替换频率的平均值与标准差高于HypoHPR区域，似乎意味着HypoHPR的突变可能导致蛋白质结构的改变，进而影响其功能^［23］。这种变化为细胞机能所不允许，部分位点的突变受到清除作用，而HyperHPR的一些位点的突变则清除压力较小，从而产生碱基替换频率的幅度波动比HypoHPR更大。这进一步说明，mtDNA的异质性受到纯化选择的影响。

通过进一步计算HyperHPR和HypoHPR中的T_S/T_V，以此作为衡量非编码区选择压力的大小。T_V对DNA双链小沟宽度和卷曲的影响比T_S大，也会影响DNA与转录因子的结合^［24］，因而，T_V更倾向于改变D-loop区原有的构象和亲和力。在完全没有选择的条件下，碱基替换完全决定于化学途径时T_S/T_V = 0.5^［25］。本研究发现，梅花鹿D-loop区HyperHPR中T_S/T_V在0～3，平均值为0.11；HypoHPR则在0～5，平均值为0.19（表3）。这说明HyperHPR更能容忍较多的T_V，而HypoHPR中的T_V受到更多的清除效应，意味着受到更强的选择压力^［26］。

由此可以看出，mtDNA的异质性是在突变和选择这对矛盾的动态之下形成的，但是也存在着显著的个性化特征。表3显示，个体平均异质性位点数在整个D-loop区的标准差达到6.77，由此可以计算得到变异系数（C_V）为0.61，而HyperHPR的C_V达到1.80，HypoHPR更是高达3.19。虽然我们目前的研究无法给出确切的原因，但也强烈提示不同个体对mtDNA功能的需求有明显的差异，这种个体差异反映在HDα上（图2），即无论是D-loop区、HyperHPR还是HypoHPR，HDα取值范围均处于较宽的置信区间，离散程度大。

异质性是由于合成mtDNA的聚合酶Polγ忠实性相对较低而导致mtDNA复制过程和损伤修复过程中产生的错误突变^［27-28］。随着年龄的增长，mtDNA的复制总量不断积累，同时因为各种因素导致的DNA损伤也不断增加^［29-30］，这一渐进性过程可使mtNDA的异质性随着年龄的增加而增加，这与本研究中检测到个体异质性频率和年龄之间显著相关的结果相一致（图3），且该现象曾在文献［31-32］中得到了证实。

综上所述，本研究揭示了梅花鹿粪便mtDNA D-loop区异质性的特点及其与年龄之间的关系。虽然样本数量还比较有限，但仍然为进一步探讨这一关系提供了重要线索。今后可以利用粪便这一采集方便的非损伤性样品，进一步挖掘来自肠道上皮细胞的mtDNA异质性信息，用于解析物种适应性演化的机制、个体对环境的适应方式^［33］、个体的适合度和动物与肠道微生物及环境（食物条件）的关系，甚至通过粪便mtDNA的异质性来估算动物的年龄。

参考文献

[1] 巫小倩，张顺华，朱砺.线粒体DNA异质性［J］.中国生物化学与分子生物学报，2017，33（1）：11-16.

[2] WU X Q， ZHANG S H， ZHU L.Mitochondrial DNA heteroplasmy［J］.Chinese Journal of Biochemistry and Molecular Biology，2017，33（1）：11-16.

[3] CALVO-RODRIGUEZ M， BACSKAI B J.Mitochondria and calcium in Alzheimer’s disease：from cell signaling to neuronal cell death［J］.Trends in Neurosciences，2021，44（2）：136-151.

[4] MANDERS F， VAN DINTER J， VAN BOXTEL R.Mutation accumulation in mtDNA of cancers resembles mutagenesis in normal stem cells［J］.iScience，2022，25（12）：105610.

[5] 赵静，万绍贵，鲍登克.线粒体损伤相关模式分子与宿主免疫调节［J］.中国生物化学与分子生物学报，2017，33（1）：38-43.

[6] ZHAO J， WAN S G， BAO D K.Mitochondrial damage-associated molecular patterns and immune regulation［J］.Chinese Journal of Biochemistry and Molecular Biology，2017，33（1）：38-43.

[7] AMORIM A， FERNANDES T， TAVEIRA N.Mitochondrial DNA in human identification：a review［J］.PeerJ，2019，7：e7314.

[8] 任晓彤，周永恒，李天峰，等.大熊猫粪便中mtDNA拷贝数与月龄的关系［J］.野生动物学报，2022，43（3）：676-684.

[9] REN X T， ZHOU Y H， LI T F，et al.The relationship between mtDNA copy number and age in giant panda feces［J］.Chinese Journal of Wildlife，2022，43（3）：676-684.

[10] NIJTMANS L G J， HENDERSON N S， ATTARDI G，et al.Impaired ATP synthase assembly associated with a mutation in the human ATP synthase subunit 6 gene［J］.The Journal of Biological Chemistry，2001，276（9）：6755-6762.

[11] MELTON T.Mitochondrial DNA heteroplasmy［J］.Forensic Science Review，2004，16（1）：1-20.

[12] CHIAL H， CRAIG J.mtDNA and mitochondrial diseases［J］.Nature Education，2008，1（1）：217.

[13] GRADY J P， PICKETT S J， NG Y S，et al.mtDNA heteroplasmy level and copy number indicate disease burden in m.3243A>G mitochondrial disease［J］.EMBO Molecular Medicine，2018，10（6）：e8262.

[14] WANG W H， LIN L J， ZHANG Q，et al.Heteroplasmy and individual mitogene pools：characteristics and potential roles in ecological studies［J］.Biology，2023，12（11）：1452.

[15] ZHANG L， YAN J B， ZHANG C，et al.Improving intestinal inflammaging to delay aging？A new perspective［J］.Mechanisms of Ageing and Development，2023，214：111841.

[16] ZENG S Y， LIU Y F， LIU J H，et al.Potential effects of Akkermansia muciniphila in aging and aging-related diseases：current evidence and perspectives［J］.Aging and Disease，2023，14（6）：2015-2027.

[17] SZCZEPANOWSKA K， TRIFUNOVIC A.Different faces of mitochondrial DNA mutators［J］.Biochimica et Biophysica Acta，2015，1847：1362-1372.

[18] GREAVES L C， NOOTEBOOM M， ELSON J L，et al.Clonal expansion of early to mid-life mitochondrial DNA point mutations drives mitochondrial dysfunction during human ageing［J］.PLoS Genetics，2014，10（9）：e1004620.

[19] WAITS L P， PAETKAU D.Noninvasive genetic sampling tools for wildlife biologists：a review of applications and recommendations for accurate data collection［J］.Journal of Wildlife Management，2005，69（4）：1419-1433.

[20] CUI L Y， LIU B Y， LI H M，et al.A simple and effective method to enrich endogenous DNA from mammalian faeces［J］.Molecular Ecology Resources，2024，24（4）：e13939.

[21] 国家重点保护野生动物名录（2021年2月1日修订）［J］.野生动物学报，2021，42（2）：605-640.

[22] List of national key protected wildlife in China （revised on February 1，2021）［J］.Chinese Journal of Wildlife，2021，42（2）：605-640.

[23] NEI M.Genetic distance between populations［J］.The American Naturalist，1972，106（949）：283-292.

[24] LIU K， XIE N， WANG Y X.Quantifying mitochondrial heteroplasmy diversity：a computational approach［J］.Molecular Ecology Resources，2024，24（1）：e13874.

[25] NEI M， LI W H.Mathematical model for studying genetic variation in terms of restriction endonucleases［J］.Proceedings of the National Academy of Sciences of the United States of America，1979，76（10）：5269-5273.

[26] 冯慧，黄原，任轶，等.陕西省林麝mtDNA D-loop区序列结构和种群遗传多样性［J］.生态学报，2014，34（20）：5887-5895.

[27] FENG H， HUANG Y， REN Y，et al.Structure of the mitochondrial DNA D-Loop region and genetic diversity of Moschus berezovskii in Shaanxi Province［J］.Acta Ecologica Sinica，2014，34（20）：5887-5895.

[28] NAUE J， HÖRER S， SÄNGER T，et al.Evidence for frequent and tissue-specific sequence heteroplasmy in human mitochondrial DNA［J］.Mitochondrion，2015，20：82-94.

[29] GUO C， MCDOWELL I C， NODZENSKI M，et al.Transversions have larger regulatory effects than transitions［J］.BMC Geno-mics，2017，18（1）：394.

[30] 任晓彤.大熊猫粪便mtDNA拷贝数和异质性特点的研究［D］.哈尔滨：东北林业大学，2022.

[31] REN X T.Characteristics of copy number and heteroplasmy of fecal mitochondrial DNA of the giant panda［D］.Harbin：Northeast Forestry University，2022.

[32] 唐晓芬，陈莉，马玉韬.密码子使用偏性量化方法研究综述［J］.基因组学与应用生物学，2013，32（5）：660-666.

[33] TANG X F， CHEN L， MA Y T.Review and prospect of the principle and methods quantifying codon usage bias［J］.Genomics and Applied Biology，2013，32（5）：660-666.

[34] BAKER A M， CERESER B， MELTON S，et al.Quantification of crypt and stem cell evolution in the normal and neoplastic human colon［J］.Cell Reports，2019，27（8）：2524.

[35] ZHENG W M， KHRAPKO K， COLLER H A，et al.Origins of human mitochondrial point mutations as DNA polymerase γ-mediated errors［J］.Mutation Research： Fundamental and Molecular Mechanisms of Mutagenesis，2006，599（1/2）：11-20.

[36] ARYAMAN J， JOHNSTON I G， JONES N S.Mitochondrial he-terogeneity［J］.Frontiers in Genetics，2019，9：718.

[37] PINTO M， MORAES C T.Mechanisms linking mtDNA damage and aging［J］.Free Radical Biology and Medicine，2015，85：250-258.

[38] HÜBNER A， WACHSMUTH M， SCHRÖDER R，et al.Sharing of heteroplasmies between human liver lobes varies across the mtDNA genome［J］.Scientific Reports，2019，9（1）：11219.

[39] DING J， SIDORE C， BUTLER T J，et al.Assessing mitochondrial DNA variation and copy number in lymphocytes of ~2，000 Sardinians using tailored sequencing analysis tools［J］.PLoS Gene-tics，2015，11（7）：e1005306.

[40] 李晓达，任晓彤，高雨薇，等.大熊猫粪便中线粒体基因库特点研究［J］.野生动物学报，2025，46（2）：269-280.

[41] LI X D， REN X T， GAO Y W，et al.Study on the characterization of the fecal mitogene pool of giant panda［J］.Chinese Journal of Wildlife，2025，46（2）：269-280.

基金资助

国家自然科学基金项目(32170517)