马淑娟,徐祎洁,何 珂,马瑞丰,朱 英*

(1.西南民族大学青藏高原研究院,成都 610041;2.浙江农林大学动物科技学院/动物医学院,杭州 311300)

Toll样受体(Toll-like receptors,TLRs)参与机体天然免疫识别、抵御外来微生物的入侵,在机体的先天免疫系统中起着重要作用[1-2]。TLR基因可分为非病毒性TLRs(non-viralTLRs)和病毒性TLRs(viralTLRs)。前者位于细胞膜上,主要识别微生物细菌的鞭毛蛋白、脂多糖、肽聚糖等[3-4];后者位于细胞器膜上,主要识别病毒未甲基化的CpGDNA和单链或双链RNA[5-6]。由于TLR直接位于宿主和环境作用的交界处,其可能与病原微生物共进化,故TLR基因是研究脊椎动物体内由病原微生物驱动的自然选择的最佳候选基因之一。

在哺乳动物中已发现13种TLRs(TLR1-13),TLR1-TLR10在陆生和水生哺乳动物中均有报道,但是TLR11-TLR13常见于小鼠、树鼩等啮齿类、树鼩类动物[6-8]。根据TLR系统发生树的相对关系,可以将TLR基因分为不同的亚家族[9-10],但哺乳动物不同类群的亚家族情况因物种数目不同而不同[11-13]。以往研究表明,所有TLR基因均受到纯化选择作用,然而不同TLR基因受到的选择压力不同[6,10],如非病毒性TLRs受到的选择压力强于病毒性TLRs。

反刍动物(偶蹄目反刍亚目)是有蹄类中物种最繁盛的动物类群,分布广、数量大,拥有200多个现生物种[14-16],包含经过人类驯化的重要家畜资源[17-20],如家牦牛、绵羊、山羊等。目前,关于反刍动物的Toll样受体基因家族的研究大部分都见于反刍动物单一物种,如对牦牛[18]、牛[21]、绵羊[21]、山羊[12]等个体的研究,或只针对TLRs中的某个或几个基因成员进行TLR遗传多样性、分子进化分析[22-24]。在反刍动物类群中,TLR家族是否存在物种特异的膨胀、TLR亚家族的分类、反刍动物不同TLR基因和不同类群反刍动物受到的选择压力模式和强度,以及反刍动物非病毒性和病毒性TLR承受的选择压力仍不清楚。

研究表明,物种数量多且覆盖全面的TLR研究能够给出较真实的系统发育关系[10]。为全面探寻反刍动物TLRs基因家族的进化模式,本研究基于公共数据库中53种反刍动物(涵盖6个科47个属)的全基因组,采用blast对反刍动物TLRs基因家族成员(TLR1-10)进行鉴定,对1)染色体位置和基因结构;2)亚家族分类;3)进化关系和选择压力;4)表达模式进行分析。研究结果有助于理解反刍动物TLR基因家族的进化,也为绵羊、牛等重要的畜牧资源和其它动物的免疫分子机制和抗病育种提供基础数据。

1 材料与方法

1.1 研究物种

反刍动物包含鼷鹿科(Tragulidae)、麝科(Moschidae)、长颈鹿科(Giraffidae)、鹿科(Cervidae)、牛科(Bovidae)和叉角羚科(Antilocapridae)6个科,76个属,拥有233多个现生物种[14-16]。在美国国家生物技术信息中心(National Center for Biotechnology Information,NCBI,https://www.ncbi.nlm.nih.gov/,Bethesda,MD,United States)上检索到117个反刍动物的基因组数据(截止到2022年12月15日),涵盖47个属。在Bos、Cervus、Muntiacus属中,共有8个物种具有三代参考基因组,这8个物种的参考基因组均保留在本研究的分析中,Moschus属仅有的2个物种(二代参考基因组)也全部保留在分析中。而在其他43个属中各选择一个组装质量高的二代参考基因组进行分析。最后,对53个物种进行了后续的分析,所选物种数目占反刍动物的22.7%(53/233),覆盖了反刍动物属的65.8%(47/76)(表1)。

表1 反刍动物TLRs基因组数据信息Table 1 Genomic data information of ruminant TLRs

1.2 TLRs基因家族成员鉴定

在NCBI上下载反刍动物TLRs的参考基因组序列,以水牛(Bubalusbubalis)的TLR1-TLR6、TLR8-TLR10(HQ327989.1、DQ288130.1、HQ343417.1、HQ343416.1、HQ327990.1、HQ327992.1、GQ499855.1、HQ242779.1、HQ327991.1),家牛(Bostaurus)的TLR7基因(EF583900.1)作为qurey序列,用BlastN进行检索,E-value为10-5,保留满足于以下两个条件的序列进行后续的分析:1)序列长度大于最长序列CDS的80%;2)不存在提前终止子。

1.3 TLRs基因家族染色体定位、基因结构分析

根据基因组数据的测序质量、完整度,选择参考基因组组装到染色体且具有注释文件的物种对其TLRs基因家族成员在染色体上的分布和基因结构进行分析。本研究对满足条件的牛科的水牛(Bovidae,Bubalusbubalis)和鹿科的马鹿(Cervidae,Cervuscanadensis),使用TBtools[25]进行TLRs基因家族在染色体上分布和基因结构的可视化。

1.4 系统发育树

在NCBI下载人类和蓝鲸的TLRs基因家族作为外类群(Human (Hosa):TLR1 NM_003263.4、TLR2 NM_001318787.2、TLR3 NM_003265.3、TLR4 NM_003266.4、TLR5 NM_003268.6、TLR6 NM_001394553.1、TLR7 NM_016562.4、TLR8 NM_016610.4、TLR9 NM_017442.4、TLR10 NM_001017388.3; Bluewhale (Blwh):TLR1 XM_036853149.1、TLR2 XM_036853661.1、TLR3 XM_036838879.1、TLR4 XM_036856405.1、TLR5 XM_036836501.1、TLR6 XM_036853144.1、TLR7 XM_036839897.1、TLR8 XM_036839825.1、TLR9 XM_036869686.1、TLR10 XM_036853063.1),使用MEGA X[26]软件对反刍动物、人类、蓝鲸的TLRs序列进行MUSCLE多重序列比对,基于最优模型(Tamura-Nei),采用最大似然法(maximum likelihood,ML)构建TLRs基因家族的系统发育树,Bootstrap值设为1 000,替换类型为核苷酸,其余参数为系统默认值[27-28]。

1.5 TLRs基因家族的选择压力分析

本研究使用非同义替换率(dN)与同义突变率(dS)的比值ω(ω=dN/dS)进行TLRs基因家族的选择压力分析。当ω>1时,表示受到正选择(positive selection);ω=1时,为中性选择(neutral selection);ω<1受到纯化选择(purifying selection)。采用MEGA X[26]计算ω。对于每个TLR基因的正选择位点,采用在线网站(datamonkey adaptive evolution server(http://www.datamonkey.org/)中的MEME(mixed effects model of evolution)、SLAC(single-likelihood ancestor counting)、FEL(fixed effects likelihood)、FUBAR(fast unconstrained bayesian approximation)4种方法进行预测。使用MEME、SLAC、FEL方法筛选正选择和纯化选择的标准是P-value<1,FUBAR方法则用后验概率>0.9进行筛选。

在上述基础上,本研究检测牛科和鹿科,病毒性和非病毒性TLR这两组的ω和正选择位点数目(positively selected sites, PSS)是否存在显著差异。利用R4.1.2的shapiro.test和var.test函数分别进行正态分布和方差齐性检验[29-30],对于符合条件的数据采用独立样本T检验,反之则采用kruskal检验[31]。

1.6 TLRs基因家族在不同组织的表达模式分析

本研究根据反刍动物基因组数据库(http://animal.nwsuaf.edu.cn/code/index.php/Ruminantia)中的RNA seq数据分析TLRs基因家族在不同组织的表达模式。反刍动物基因组数据库包含多个物种的RNA seq数据,由于绵羊(Ovisaries)具有丰富的组织表达数据,故本研究采用绵羊的RNA seq进行后续分析。使用TBtools软件绘制其内分泌系统、免疫系统、神经系统、呼吸系统、消化系统共25个组织的表达模式热图。

2 结 果

2.1 反刍动物TLRs基因家族的数量

对53个反刍动物的全基因组进行TLR鉴定,得到53种反刍动物(鼷鹿科 1个、叉角羚科 1个、长颈鹿科2个、鹿科 11个、麝科 2个、牛科 36个,表1和图1)共530条TLR基因序列。对于每一个TLR基因都仅鉴定出1条TLR序列,不存在多拷贝现象。通过去除短序列(67条)、假基因序列(42条),剩余421条TLR真基因序列,TLR1-TLR10的10个基因分别有39、40、45、43、44、36、43、41、44、46条序列(表1)。在科水平上,鼷鹿科共鉴定到8个TLR基因(无TLR3、TLR10),叉角羚科鉴定到9个TLR基因(无TLR8),长颈鹿科、鹿科、麝科、牛科都鉴定到10个TLR基因(表1和图1)。在物种水平上,通过剔除短序列和假基因之后,鉴定到的TLR基因数目为0～10个。牛科的高角羚(Aepycerosmelampus)和柯氏犬羚(Madoquakirkii)鉴定到的TLR真基因数目为零。

进化树改自文献[14]的图1Phylogenetic tree was modified from figure 1 in reference [14]图1 反刍动物的TLR基因数目Fig.1 TLR genes number of ruminants

2.2 反刍动物TLRs的染色体定位、基因结构

水牛和马鹿的TLRs家族成员在染色体上的定位不同、基因结构不同(图2)。水牛TLR1、TLR6、TLR10连锁,并位于7号染色体,TLR7和TLR8连锁,位于X染色体,TLR2、TLR3、TLR4、TLR5、TLR9分别分布于17、1、3、5、21号染色体上(图2a)。马鹿TLR1、TLR6、TLR10位于19号染色体,TLR7和TLR8位于X染色体,TLR2、TLR3、TLR4、TLR5、TLR9分别分布于染色体1、31、30、13、22上(图2b)。

a.水牛TLRs基因在染色体上的定位;b. 马鹿TLRs基因在染色体上的定位;c. 水牛TLRs基因结构;d. 马鹿TLRs基因结构a. Chromosomal localization of TLR genes in Bubalus bubalis; b. Chromosomal localization of TLR genes in Cervus canadensis; c. TLR genes structure of Bubalus bubalis; d. TLR genes structure of Cervus canadensis图2 水牛和马鹿TLRs基因在染色体上的定位和基因结构Fig.2 Chromosomal localization and gene structure of TLRs in Bubalus bubalis and Cervus canadensis

水牛和马鹿的CDS全长是保守的,但基因结构分析表明,水牛和马鹿在外显子数目、基因长度存在差异(图2c,图2d)。水牛的TLR基因含有2～5个外显子,其中TLR3基因有最多的外显子数目,水牛的TLR基因长度位于7 828 (TLR9) ～33 883 bp (TLR10)之间(图2c)。马鹿TLR基因含有2～7个外显子,其中TLR7基因具有最多的外显子。马鹿TLR基因基因长度在5 119(TLR9)～136 030 bp (TLR2)之间(图2d)。不同基因的长度差异来自于内含子长度的差异上。

2.3 反刍动物的TLRs的分类

以蓝鲸和人类TLRs基因作为外类群,与421条反刍动物TLRs基因序列构建的ML系统发育树(图3)表明,反刍动物TLR序列聚成5个基因簇,分别是TLR1 (TLR1/6/10)、TLR2、TLR4、TLR3(TLR3/5)、TLR7 (TLR7/8/9)亚家族。TLR1和TLR7亚家族拥有最多的TLR家族成员,TLR2和TLR4仅有1个直系同源基因。在每个TLR基因家族成员内部,TLR基因序列按照科聚类,叉角羚科和长颈鹿科、麝科和牛科各自形成姐妹群(图3)。

自展值大于70%的用*表示Bootstrap values larger than 70% are denoted by * 图3 反刍动物TLR基因的最大似然树Fig.3 ML tree of ruminant TLR genes

2.4 反刍动物TLRs自然选择压力

本研究检测了53种动物10个TLR基因的自然选择压力模式和强度,非同义替换率(dN)和同义替换率(dS)的比值(ω)表明所有的TLR基因均受到纯化选择作用(表2),ω值从TLR8的0.121到TLR6的0.369,提示反刍动物的TLR基因受到一定的功能约束。

表2 反刍动物TLRs基因家族进化压力Table 2 Evolutionary pressure of TLRs gene family in ruminants

尽管反刍动物的TLR基因整体受到明显的纯化选择作用,但是在除了TLR7外的其他的TLR基因中均检测到了显著的正选择位点(193个),其中35个是4种方法均鉴定到的(表2),提示这些位点受到强烈的正选择作用。不同TLR基因的正选择位点数目差异明显,从TLR7的0到TLR1的6.80‰。在所有的直系同源TLR基因中,TLR8和TLR7分别拥有最低的ω值和最低比例的正选择位点,提示这两个基因受到的选择约束高于其他直系同源TLR基因。牛科和鹿科ω均小于1,都受到强烈的纯化选择,但是相比于鹿科,牛科具有更高的ω值和更高比例的正选择位点,提示鹿科TLR基因受到的约束力高于牛科(图4a-b,表2;牛科ω=0.253,鹿科ω=0.218, χ2=7.303, df =1,P=0.007;牛科PSS比例为0.253‰,鹿科PSS比例为0.218‰,χ2=0, df=1,P>0.05)。在所有反刍动物的TLR基因中,非病毒性TLRs(TLR1/2/4/5/6/10)比病毒性TLRs(TLR3/7/8/9)具有更高的ω和更多的正选择位点数目(图4c-d,病毒性ω=0.188,非病毒性ω=0.264, t=1.650, df=8,P=0.138;病毒性PSS比例为2.67‰,非病毒性PSS比例为5.20‰, t=2.586, df=8,P=0.032)。

a. 牛科和鹿科的ω; b. 牛科和鹿科的正选择位点的比例; c. 非病毒性和病毒性TLR的ω; d. 非病毒性和病毒性TLR的正选择位点的比例。ns. 差异不显著a. ω of Bovidae and Cervidae; b. Percentage of positively selected sites in Bovidae and Cervidae; c. ω of non-viral and viral TLRs; d. Percentage of positively selected sites in non-viral and viral TLRs. ns. Not significant图4 反刍动物TLR基因的选择压力Fig.4 Selection pressure of TLR genes in ruminants

2.5 绵羊TLRs基因家族在不同组织的表达模式

本研究对绵羊10个TLRs基因家族成员在不同系统和组织中的组织表达量计算算术平均值,通过表达量热图发现几乎所有TLR基因都在免疫组织器官上表达(图5,灰色色块颜色所示)。TLR2/9/1/4在免疫组织的表达量整体高于TLR3/7/5/8/6/10,TLR2和TLR9基因在外周血单个核细胞(peripheral blood mononuclear cell,PBMC)中存在高水平的表达。

图5 绵羊TLR基因表达模式Fig.5 Expression patterns of TLR genes in Ovis aries

3 讨 论

反刍动物对于人类社会具有重要的价值及生物学意义,对于反刍动物的TLR家族研究有助于理解反刍动物免疫系统的进化机制。在已报道文献中,通过PCR扩增的方法确定牦牛、山羊、牛、绵羊[12,18-19,32]具有TLR1-TLR10 10个TLR基因,这些研究仅包含了13种反刍动物,约占反刍动物总数的5.6%(13/233)。本研究借助于NCBI中丰富的反刍亚目基因组数据,通过BLAST对反刍亚目的6个科,47个属(约占65.8%),53个物种(约占物种总数的22.7%)的TLR基因家族进行了鉴定。较以往研究,本研究的物种覆盖了反刍亚目全部6个科及一半以上的属,首次对反刍亚目谱系的TLR基因家族的基因拷贝数目、自然选择压力进行了研究。

本研究发现,反刍亚目的10个TLRs基因均是单拷贝,表明TLR在反刍亚目中相对保守。在其他反刍动物的研究中,如牦牛[18]、山羊[19]、水牛[12],TLR1-10基因均为单拷贝,与本研究结果一致。在鸟类研究中发现TLR1、TLR2存在多拷贝现象[3]。除了鸟类,脊椎动物的其他类群都具有TLR1-10基因,表明TLR基因数目在脊椎动物中也是保守的。本研究采用小鼠的TLR11/12/13基因作为query序列,在反刍亚目所有物种中鉴定到的序列长度过短(低于query序列的50%),推测TLR11/12/13在反刍亚目中可能为假基因。对于TLR1-10,某些物种鉴定到的假基因,根据近源物种的情况,推测是由测序或者组装拼接错误导致的。

本研究基于421条TLR序列的系统发育树,将反刍亚目的TLR基因分为5个亚家族,分别是TLR1家族(TLR1/6/10)、TLR2家族(TLR2)、TLR3家族(TLR3/5)、TLR4家族(TLR4);TLR7家族(TLR7/8/9)。这与以往的TLR基因进化研究结果不完全一致[6,10,12]。在已有研究中,TLR1/6/10由于具有较高的序列相似性都被归为同一个亚家族,TLR7/8/9与此类似,也聚为一簇,但TLR4与其他9个TLRs基因差异较大,单独聚为一簇,与本研究结果一致[9-10]。

TLR1亚家族是最大的亚家族之一,包含TLR1、TLR6、TLR10。水牛、马鹿(本研究)、家牛[33]、人类[6,10]的染色体位置表明TLR1、TLR6、TLR10位于同一条染色体,这一结果支持了3个TLRs基因进化可能是非独立的。此外,这3个基因的高度相似性可能源于配体结合的特性,如TLR1、TLR6、TLR10与TLR配体结合作为异二聚体以传递细胞信号并诱导相关的免疫反应[34]。而TLR7亚家族的TLR7/8/9都在细胞内膜上表达,都具有种属特异的配体,都属于病毒型TLRs基因。有研究表明,TLR8和TLR9与它们的配体结合作为内体[35],表明保守的TLR在组合和传输细胞信号的过程中通过相似的模式在组合配体中起重要作用。此外,TLR7和TLR9在病毒衍生基序的TLR中的选择模式与细菌PAMPs相互作用的模式有很大不同。而TLR7和TLR8位于同一条染色体。其他未组装到染色体的基因组数据表明,超过80%的物种都显示TLR1、TLR6和TLR10以及TLR7和TLR8都分别在同一个scaffold。以上可能是TLR1/6/10和TLR7/8/9分别被归为同一亚家族的原因。TLR4亚家族只有TLR4基因,该基因是识别革兰氏阴性菌LPS的一个免疫基因[34],与其他TLR基因相似性较低,表明TLR4与其他TLR基因可能在功能上存在差异。

而TLR2、3、5在以往不同的研究中聚类情况不尽相同。TLR2在有些研究中与TLR1/6/10被纳入TLR1亚家族,如Liu等[10]对脊椎动物TLRs基因家族的研究。在本研究中,TLR2与TLR1/6/10聚为不同的基因簇,提示在反刍动物的进化过程中TLR2基因与其他3个基因产生了较大的分化。在脊椎动物TLR进化中,TLR3和TLR5聚为不同的簇,而反刍动物的TLR3和TLR5在序列上更为相似,故聚为一支,提示TLR3和TLR5在反刍动物外产生了分化。

ω结果表明,反刍动物TLRs基因家族成员都经历了纯化选择,这与两栖类[5-6]、鸟类[3,36]、鱼类[37]等其它脊椎动物研究结果一致,表明TLRs基因功能在进化上是相对保守的。值得一提的是,反刍动物的TLRω(0.121～0.369)与快速进化的基因,如极光激酶(AKs)[38]、DAZ基因类似[39]。在科水平上,牛科总体选择压力大于鹿科,正选择位点数量也多于鹿科反刍动物,表明鹿科在进化水平上受到更强的约束力。病毒型TLRs比非病毒型TLRs受到更强的功能选择约束,这与前人关于灵长类、树鼩类、两栖类动物研究的结果一致[6,13]。这可能因为病毒TLRs在免疫应答中维持其识别病毒核酸功能,同时也要对自身损伤的结构进行识别,从而避免非同义替换的出现或者非同义突变的积累,以此来维持自身的功能[34]。此外,非病毒TLRs有些在功能上会存在冗余,会识别相同的病原微生物PAMPs,能够接受该类病毒发生非同义突变,在种群受到正选择时保留下来[6,34]。

绵羊的TLRs基因在免疫系统组织中高表达,该结果与以往研究结果一致[7,12,18],表明TLRs在抵抗病原体入侵中发挥重要作用。TLR2和TLR9基因在PBMC组织上高表达,PBMC包含较多的淋巴细胞(T细胞、B细胞、NK细胞),是先天免疫系统和适应性免疫系统的关键组成部分,可以很好的保护机体免受病毒、细菌的感染[40]。TLR3基因在绵羊内分泌系统中存在高表达现象,这可能与其具有独特的识别双链RNA病毒的功能有关,表明TLR3在绵羊内分泌系统疾病的发生中发挥重要作用,可能为绵羊生殖疾病引起的生殖能力低下和其它问题的解决提供有效帮助[7,34]。此外,本研究发现绵羊的TLR1和TLR4表达相近,TLR6和TLR10在所有组织中表达水平相似,这可能是因为它们在鉴定病原微生物和参与生物体的抗病性和免疫应答过程中相互配位且共同发挥作用[7-8]。

4 结 论

本研究阐明了反刍动物TLR基因家族的基因拷贝数目、系统发育关系、选择压力。结果表明,反刍动物的TLR1-10均是单拷贝基因,可归为5个亚家族。尽管每个TLR基因在整体上都受到强烈的纯化选择作用,但是仍有氨基酸位点受到正选择作用。相对于牛科,鹿科在进化历程中受到更强的选择约束。反刍动物的病毒性TLRs比非病毒性TLRs受到更强的自然选择约束。由于本研究的局限性,无法确认在某些物种中鉴定到的TLR1-10假基因确实是由测序造成的。由于鼷鹿科、长颈鹿科和麝科所包含的物种数目均少于3个,不能满足正选择位点的分析条件,故这3个科在进化过程中的所受到的选择作用仍不清楚。此外,虽然本研究利用多种方法检测到反刍动物的多个正选择位点,但是缺少疾病的表型数据,这些正选择位点是否与疾病相关联仍然是未知的,从而不能更好地反映Toll样受体蛋白与其致病对应物之间的进化关系。在未来的研究中,可以通过正选择位点与疾病的关联分析并结合正选择位点在蛋白三维结构功能区的定位研究来确认正选择位点的功能。