不同品系豚鼠MHCⅠ类基因的多态性及差异表达

2021-09-10卫振何珂洪胜辉刘迪文

浙江大学学报（农业与生命科学版） 2021年4期

卫振，何珂，洪胜辉，刘迪文*

（1.浙江大学实验动物中心，杭州 310058；2.浙江农林大学动物科技学院和动物医学院，杭州 311300）

豚鼠实验动物模型，在模拟人类病理学和感染的研究中被广泛使用[1]。该物种已被用于研究心力衰竭、骨关节炎等复杂的病理过程。其中，豚鼠被用于口蹄疫病毒（foot-and-mouth disease, FMD）研究的时期长达数十年[2]。FMD 是一种人畜共患病，严重危害畜牧业和人类健康。虽然豚鼠是测试FMD疫苗效力的理想模型，但由于其不同品系的遗传结构及毛色性状对疫苗检测存在影响，常常导致豚鼠对口蹄疫病毒攻毒实验呈现不同的结果，这可能使得人们难以正确评估FMD疫苗的免疫效果[3]。

2008 年，我们培育出Zmu-1∶DHP 远交系豚鼠（体白色）和Zmu-2∶DHP 远交系豚鼠（体黑色），发现两者与其亲本DHP 远交系豚鼠（体花色）在表型和遗传结构上均存在显著差异，主要表现为Zmu-1∶DHP对FMD的敏感性为100%，显著高于DHP豚鼠（亲本，0～90%的敏感性）[3]。进而用组胺激发试验检测其气道反应性，结果显示，Zmu-1∶DHP豚鼠气道平滑肌对组胺的敏感性大于DHP豚鼠[4]。而在遗传结构上，用微卫星标记分析3 个品系的遗传多样性，结果发现这3 个品系之间存在显著遗传变异[5]。由于不同品系豚鼠对病毒敏感程度存在差异，尤其是Zmu-1∶DHP品系对病毒侵染敏感，因此，有必要对豚鼠的免疫分子遗传差异进行研究。主要组织相容性复合体（major histocompatibility complex,MHC）是启动免疫防御的主要因子，也是脊椎动物获得性免疫过程中的重要分子[6]，它专门用于向T淋巴细胞提供外来抗原。早期血清学研究表明，豚鼠MHC 基因存在3 个具有差异但是连锁的基因区域[7]。尽管豚鼠被用作传染性动物模型来研究免疫遗传机制，但并未见有关其MHC基因的深入报道。在本研究中，我们解析了豚鼠MHCⅠ类基因的结构，分析了不同豚鼠品系间MHCⅠ类基因单倍型表达是否存在差异，旨在为不同品系豚鼠抗病分子标记筛选及疫苗的制备奠定理论基础。

1 材料与方法

1.1 材料

1.1.1 实验动物

Zmu-1∶DHP 远交系豚鼠（体白色，标记为W，14 只）、Zmu-2∶DHP 远交系豚鼠（体黑色，标记为B，10只）和亲本DHP品系豚鼠（体花色，标记为M，12 只），各品系豚鼠雌雄各半，均为3 月龄（体质量300～450 g）。实验动物由浙江大学实验动物中心[SCXK（浙江）2012-0052 和SCXK（浙江）2012-0178]提供，本实验得到浙江大学实验动物中心实验动物使用和管理委员会的批准（ZJU20170320）。

1.1.2 主要试剂与仪器

Trizol 试剂（美国Invitrogen 公司）、反转录酶（日本Toyobo公司）、DNA 消化酶（北京康为世纪生物科技有限公司）、Nanodrop 2000 超微量分光光度计（美国赛默飞世尔科技有限公司）、Hiseq PE300高通量测序仪（杭州联川生物技术股份有限公司）。

1.1.3 样品采集及目标区域的扩增测序

仔细解剖得到豚鼠的脾脏，在液氮中快速冷冻。用Trizol试剂从100 mg样品中提取总RNA，通过DNA 消化酶去除基因组DNA，并使用反转录酶合成cDNA。参考人、大鼠和小鼠的MHC序列设计引物（上游引物，5′-ATGGCGGAGGTGGCCCCCC GAACCCT-3′；下游引物，5′-TCTGCTCYTCTCC AGAAGGCACCA-3′）来扩增目标区域，长度为672 bp，包括部分外显子1、整个外显子2和3以及部分外显子4，其中外显子2、3含MHC Ⅰ类基因的抗原呈递区域，为高度可变区。操作如下：随机选择6个个体，以cDNA 为模板进行聚合酶链式反应（polymerase chain reaction, PCR）扩增，获得目的条带产物后将其连接到pMD18T 载体上；挑选阳性克隆（每个个体20 个）进行测序并将相关序列作为参考设计高通量测序的引物。

由于测序片段的限制，在高通量测序中利用巢式PCR进行扩增。第1轮扩增获得一段长度为460 bp的片段（Hi-up，5′-CCGAACCCTCCTCCTGCT GCTGTC-3′；Hi-dn，5′-CCGCCACCCAGGWGCT CAGGTCC-3′）；在第2轮扩增中采用分段扩增，将其分成A和B片段（A片段引物为Hi-up，5′-CCGA ACCCTCCTCCTGCTGCTGTC-3′和Hi-A-dn，5′-CCGCCACCCAGGWGCTCAGGTCC-3′；B片段引物为Hi-B-up，5′-CCGAGYGRRCCTVGGAMYCT GCKC-3′和Hi-B-dn，5′-CAGGGCCCAGCACCG CAGGGTGAC-3′）。A和B 2部分具有184 bp的重叠区（图1），能够有效进行数据拼接。待扩增完成后，将20 个个体的PCR 产物混合为一个库，进行高通量测定（使用个体标签序列来进行个体归属）。

图1 本研究所用的PCR引物示意图Fig.1 Location of PCR primers used in this study

1.2 数据分析

1.2.1 测序数据分析

按以下步骤分析高通量测序数据：1）使用个体标签序列将数据分组到对应个体中；2）选择前100个测序片段（reads）的序列用于后续的研究，并将相同序列的测序片段加在一起（如果引物部分中存在一些不同的核苷酸，将它们作为相同的等位基因）；3）使用MEGA 5软件[8]比对片段A和片段B的序列，并使用184 bp的重叠区域进行组合；4）根据上一步中的序列，从克隆测序结果中鉴定并确定单倍型，单倍型命名为Capo-ⅠA*。

1.2.2 序列进化分析

将单倍型序列翻译成氨基酸（amino acid,aa）序列后，参照HLAⅠ类基因中抗原结合位点（peptide binding site,PBS）氨基酸的定位，推断豚鼠MHCⅠ类基因的抗原结合位点。利用SWISS-MODEL（https://swissmodel.expasy.org/）预测蛋白质的三级结构。利用MEGA X 计算抗原结合域（peptide binding region, PBR）α1 和α2 的核苷酸距离和氨基酸距离。通过多种方法来检测豚鼠MHCⅠ类基因的选择作用。首先，针对整个抗原结合域、PBS及非抗原结合位点（non-PBS），选择Jukes-Cantor 修正的Nei&Gojobori 距离模型计算非同义替换率（dN）和同义替换率（dS），通过dN/dS比值来判断选择发挥的作用。其次，使用混合效应模型（mixed effects model of evolution,MEME）来计算正向选择的位点。此外，利用重组检测遗传算法（genetic algorithm for recombination detection,GARD）检查重组[9]。GARD的检测在http://www.datamonkey.org/网站上进行。

1.2.3 系统进化分析

使用邻接法（neighbor-joining method,NJ）对鉴定的单倍型与其他物种构建系统进化树。

1.2.4 品系单倍型差异分析

利用测序片段数目的不同，在不同品系间进行单倍型表达情况的比较分析。PCR 的偏好性可能会导致测序片段数目存在偏差，但在不同个体中扩增偏好性会同时存在，因此，该分析比较是可以进行且可信的。使用t检验分析不同品系的单倍型频率之间的差异显著性。

2 结果与分析

2.1 豚鼠MHCⅠ类基因的序列变异结果

通过对672 bp片段进行PCR克隆和测序，获得多条豚鼠MHCⅠ类序列，据此设计引物进行MHCⅠ类抗原结合域（α1 和α2）的高通量测序。结果显示，大多数个体的测序片段数目超过3 000，符合本研究的分析要求，另外有3 个个体的测序片段数目低于3 000（附表1，http://www.zjujournals.com/agr/CN/10.3785/j.issn.1008-9209.2020.10.222），因此，未加入后续分析。共使用33个个体进行数据分析，平均测序片段数目为10 731，足以说明等位基因信息和表达情况。

通过比较重叠区域184 bp 片段中的独特核苷酸，发现片段A 中的所有等位基因都可以和片段B完全匹配，并且所有拼接好的序列（片段A＋片段B）均能在克隆测序中得到确定（附表2，http://www.zjujournals.com/agr/CN/10.3785/j.issn.1008-9209.2020.10.222）。因此，我们采用单倍型数据（包括α1 和α2 区域，GenBank 登录号分别为MH841957～MH841962 和MH841964～MH841965）进行下一步分析，未加入分析的个体中未能找到新的单倍型序列。单倍型片段长度为673或603 bp，可编码224或201 个氨基酸，包括部分的外显子1（18 个氨基酸）、完整的外显子2（α1区域，90个氨基酸）、完整的外显子3（α2 区域，92 或69 个氨基酸）和部分的外显子4（24个氨基酸）。整个DNA序列中没有发现提前终止密码子，但单倍型Capo-ⅠA*06 存在69 bp 的片段缺失（图2）。

从图2 中可知，8 种单倍型中存在88 个氨基酸多态性位点，其中α1 区域和α2 区域分别有43 和45个多态性位点。单倍型Capo-ⅠA*04 和Capo-ⅠA*05的α2区域具有相同的氨基酸序列。PBS和抗原呈递区域存在密切关系。研究结果发现，豚鼠MHCⅠ类基因单倍型的58个PBS中存在41个多态性位点，且Capo-ⅠA*06和Capo-ⅠA*07具有相同的PBS（不包含Capo-ⅠA*06缺失的位点）。因此，在后期分析中可考虑将这些单倍型进行合并处理。

图2 豚鼠MHCⅠ类基因单倍型多态性分析Fig.2 Analysis of haplotype polymorphism of MHC classⅠgene in guinea pig

通过对蛋白结构进行预测，发现该缺失片段位于α2 区域的起始处，该区域含有关键的第11 位胱氨酸，这个位点与二硫键的形成有关（图3）。

图3 单倍型Capo-ⅠA＊01和Capo-ⅠA＊06的蛋白质三维结构预测图Fig.3 Protein three-dimensional predicted structures of Capo-ⅠA*01 and Capo-ⅠA*06

2.2 豚鼠MHCⅠ类基因单倍型的进化选择

非同义替换率（dN）表示核苷酸突变导致相应位点的氨基酸发生了改变，而同义替换率（dS）表示不产生氨基酸的变化，二者的比值（dN/dS）可用来衡量功能基因是否受到选择作用。通常认为，当dN/dS＞1时，该基因受到正向选择作用，而当dN/dS＜1时，该基因受到净化选择作用。我们分别计算了豚鼠MHCⅠ类基因单倍型α1 和α2 的dN/dS值。结果发现，在α1和α2这2个区域均存在明显的正向选择作用（dN/dS值分别为2.108 和2.509），利用整个单倍型计算dN/dS存在类似的结果（dN/dS=2.319）。采用MEME检测得到17个正向选择作用的氨基酸位点，这些位点绝大多数（71.59%）都位于PBS 中，只有5个（29.41%）不在PBS中。

进化树分析结果（图4）表明，豚鼠的所有22 条序列[其中14 条序列来源于NCBI 数据库（https://www.ncbi.nlm.nih.gov/）]形成一个大簇。这表明本研究的8条序列能够代表数据库中豚鼠其他单倍型的进化特征，其中，Capo-ⅠA*04 和Capo-ⅠA*05聚为一支（支持率为98%），而剩下的序列聚类支持率均不高。因此，我们仅采用本研究中8条单倍型序列用于后续分析。

图4 豚鼠MHCⅠ类序列系统发育树Fig.4 Phylogenetic tree of MHC classⅠsequences in guinea pig

在加入其他哺乳动物的相关序列后，进化树聚类分析结果（图5）表明，序列基本按照物种聚类，而并非按照位点聚类。从中可知：豚鼠的8 个单倍型序列分为2大簇，分别为Capo-ⅠA*04/*05（支持率为100%）和其他绝大多数单倍型（支持率为87%），而这2支的聚类支持率为71%。此外，大鼠、人和猪的相关序列基本按照物种成簇分布，这说明在哺乳动物的MHCⅠ类基因中存在协同进化的现象。

图5 部分哺乳物种MHCⅠ类序列的系统发育树Fig.5 Phylogenetic tree of MHC classⅠsequences in some mammal species

通过序列对比发现，单倍型之间很多区段对应的核苷酸序列几乎完全一致，因此，推断在豚鼠的MHCⅠ类多个位点间可能存在密切的基因交换或重组事件。利用GARD方法检测Capo-ⅠA单倍型的重组，结果显示，外显子2（α1）的第71位氨基酸是一个重要的重组断点（P＜0.001）。虽然在α2 中未检测到重组事件，但多个单倍型之间共享DNA 片段，如Capo-ⅠA*04和Capo-ⅠA*05、Capo-ⅠA*06和Capo-ⅠA07，表明豚鼠的单倍型多样性可能在基因重复、基因重组或协同进化的作用下形成。

2.3 不同品系间豚鼠MHCⅠ类单倍型表达差异

在3 个不同品系中检测了MHCⅠ类单倍型的表达情况。结果（表1）发现，在每个品系中都存在这8条单倍型（外显子2上存在8条等位基因），因此，推测在3个品系个体中都至少存在4个MHCⅠ类表达位点。根据3 个品系中单倍型的表达情况，分析发现：Zmu-2∶DHP 和DHP 品系仅Capo-ⅠA*06 测序片段频率存在显著差异（P＜0.01）；Zmu-1∶DHP 和DHP品系存在3个单倍型表达差异，分别是Capo-ⅠA*01（P＜0.01）、Capo-ⅠA*04（P＜0.01）和Capo-ⅠA*05（P＜0.01），以及Zmu-1∶DHP 和Zmu-2∶DHP品系存在4个单倍型表达差异，分别是Capo-ⅠA*01（P＜0.01）、Capo-ⅠA*03（P＜0.05）、Capo-ⅠA*05（P＜0.01）和Capo-ⅠA*07（P＜0.05）（表1）。其中，Zmu-1∶DHP 品系和其余2 个品系之间的Capo-ⅠA*01和Capo-ⅠA*05表达频率均存在显著差异。

表1 3个品系间MHCⅠ类基因单倍型频率分析及比较Table 1 Analysis and comparison of frequencies of MHC classⅠgene haplotypes in three strains

根据之前对PBS的分析结果，将抗原呈递区域相同的单倍型Capo-ⅠA*06 和Capo-ⅠA*07 的测序片段合并在一起进行显著性比较，发现Zmu-1∶DHP 和DHP、Zmu-2∶DHP 品系之间均差异显著（P＜0.05），而DHP与Zmu-2∶DHP品系之间不存在显著性差异。

3 讨论

本研究利用高通量混合池文库测序的方法，首次系统地研究了不同品系豚鼠的MHCⅠ类基因多态性，分析了α1和α2区域的选择作用、重组和系统进化，并且调查了不同品系豚鼠中单倍型的表达情况，结果如下：1）在调查的豚鼠个体中，至少存在8条表达的MHCⅠ类基因单倍型，因此，推测至少存在4个MHCⅠ类基因表达位点；其中单倍型Capo-ⅠA*06在外显子3（α2区域）缺少23个氨基酸；8条单倍型序列中多态性位点达到88个；2）通过dN/dS及MEME检测，发现豚鼠MHCⅠ类基因抗原呈递区域受到显著正向选择作用，而非抗原呈递区域受到净化选择；3）在不同品系间进行单倍型表达情况分析，发现Zmu-1∶DHP品系的单倍型表达存在显著变化，和其他2个品系相比，存在多个表达单倍型测序片段频率差异，推测这可能和品系的差异免疫能力有关。

3.1 MHC 进化

MHC作为和免疫密切相关的基因家族，在进化速度上较其他基因快[10]，具有高度多态性。在长期的物种进化过程中，MHC 基因家族经历了各种变化，如基因重组、染色体分裂、基因重复等[11]，才形成现在适合各个物种的既保守又特异的结构。目前，研究认为HLAⅠ类基因有6 个经典位点和10 个非经典位点，而在其他哺乳动物中也同样存在MHCⅠ类基因的经典位点和非经典位点[12]。非经典位点与经典MHCⅠ类基因的差异在于多态性的高低和其编码产物的组织分布广泛度。在分析过程中我们发现，豚鼠个体中最多存在8条表达的MHCⅠ类单倍型，及单个外显子上存在8条等位基因，因此推测至少存在4 个MHCⅠ类表达位点。但由于未能通过内含子等信息来判断定位位点信息，系统进化树也无法通过单倍型之间的差异推断出相关的位点分类（仅Capo-ⅠA*04和Capo-ⅠA*05形成一支明显的进化簇），因此，无法将这8 个单倍型分配到相对应的位点。此外，多物种的MHCⅠ类序列进化树清晰地表现出该基因家族的协同进化现象。由于MHC基因家族存在重复单元，因此在位点间存在基因转化事件的发生，这会在物种间形成多个拷贝且相似性很高的多个基因，是MHC 基因家族平衡选择的重要机制[13]。

在分析过程中，我们发现Capo-ⅠA*06 单倍型的α2 区段出现23 个氨基酸的缺失，将影响α2 区域第11 位和第74 位二硫键的形成，该结构已被报道和抗原呈递有一定关系[14]。如牦牛DQA2基因α1、α2 功能区的单核苷酸多态性（single nucleotide polymorphism,SNP）改变可能会对机体细胞的抗原呈递和跨膜运输产生影响，进而影响生物体对疾病的抗性[15]。通过数据库分析发现：在多个物种[如毛丝鼠（XP_013365245）和鼹鼠（XP_01906609）]中也存在缺失部分α2区段现象；在海蟾蜍[16]、树鼩[17]和大鲵[18]上报道存在MHCⅠ类基因的可变性剪切。由此，我们推测该单倍型可能为低表达（在测序数据中得到验证，测序片段数目不超过1 000）。豚鼠的MHCⅠ类单倍型PBS 存在正向选择作用（dN/dS＞1），而在非PBS区存在净化选择作用（dN/dS＜1）。本研究采用的豚鼠均为人工品系，说明在人工选择过程中并没有弱化自然选择的作用，或者自然选择的作用得到了加强。

3.2 等位基因及单倍型差异表达

影响基因表达差异的因素有很多，主要有基因多态性和功能基因的差异表达2方面。前者由等位基因之间的核苷酸变异导致形成，后者可能由于非编码区的变异对等位基因的表达形成影响。现有研究已明确，不同等位基因编码的HLAⅠ类分子在调节免疫效应细胞功能中常常发挥不同作用。某些等位基因的产物和自然杀伤细胞上的受体分子结合，参与调节NK细胞的活化[19]；特定的等位基因产物结合特定抗原肽，与T细胞抗原受体结合，在特异性T细胞的活化中发挥关键作用[20]。有关白血病细胞的研究发现，不同等位基因编码的HLAⅠ类分子表达水平变化不同[21]。因此，在理解机体在生理和病理条件下的免疫调节机制中，需要将同一基因座位不同等位基因的表达水平区别对待，若等位基因位于不同的位点上，其发挥的作用也不同[22]。

目前，在对HLA 分析的过程中，通过等位基因特异性系列引物，已成功建立了实时定量PCR 方法，以检测同一个体2条不同的HLA-C等位基因的表达水平[23]。在本研究中，我们发现8 条单倍型在不同品系间的测序片段（reads）频率存在差异，其中Zmu-l∶DHP 品系的多个单倍型测序片段（reads）频率均显著不同于Zmu-2∶DHP 品系和DHP 品系（Capo-ⅠA*01、Capo-ⅠA*03、Capo-ⅠA*04、Capo-ⅠA*05 和Capo-ⅠA*07），而且除Capo-ⅠA*03外，其余4个单倍型测序片段（reads）数目均超过1 000，具有较高的可信度。同时，结合系统进化树分析结果，推测上述单倍型表达差异可能会造成品系间抗病能力的差异。

我们在前期的研究中发现，Zmu-l∶DHP豚鼠较适于制作FMD 病毒感染模型，DHP 豚鼠适于研究T、B淋巴细胞及生物制品的免疫功能。O-A型口蹄疫疫苗对2个品系豚鼠诱导出的特异性IgG抗体在Zmu-l∶DHP品系豚鼠中维持时间较长，而DHP豚鼠T 淋巴细胞增殖率和细胞因子分泌量显著高于Zmu-l∶DHP 豚鼠[3]。这些免疫性能上的差异，可能和豚鼠MHCⅠ类基因的单倍型表达差异有关。