从头组装豚鹿血液转录组揭示雌雄基因表达差异

刘运，柳方圆，余建秋，牛李丽，刘选珍，李静*

（1. 四川大学生命科学学院，生物资源与生态环境教育部重点实验室，成都610065；2. 成都动物园，成都野生动物研究所，成都610081）

豚鹿Axis porcinus属偶蹄目Artiodactyla 鹿科Cervidae 花鹿属Axis，包括 3 个亚种：印度亚种A. porcinus annamiticus、南亚亚种A. porcinus kuhli和指名亚种A. porcinus porcinus，其中，印度亚种分布在巴基斯坦、尼泊尔、印度、孟加拉国和缅甸，南亚亚种分布在印度，指名亚种分布在中国、泰国、老挝、柬埔寨和越南（Tanushree & Mathur，2000）。豚鹿在我国仅分布于云南省西南部，曾一度被认为已在国内灭绝（胡婧等，2011）。近年来，野生豚鹿种群数量急剧减少，自2008 年起已被世界自然保护联盟（IUCN）列为濒危（EN）物种（Timminset al.，2015）。雌雄豚鹿在外形上区别较大：雄性更大且有鹿角。研究雌雄豚鹿基因表达的特征并揭示鹿角生长发育相关基因对于理解其生物学特性具有重要意义。

通过构建系统发育树，Abbas 等（2017）发现豚鹿具有遗传多样性较低的特点；Wang 等（2017）采用重测序技术，以单核苷酸多态性（SNPs）为标记研究了成都动物园圈养豚鹿的遗传多样性也支持遗传多样性低的结论；Wang 等（2019）首次报道了豚鹿的全基因组序列，为鹿科动物比较基因组学和遗传育种等研究提供了参考。但该基因组并未对豚鹿的基因进行功能注释，豚鹿基因转录表达方面的研究尚未见报道。

转录组可对不同组织或样品间的基因表达量进行定量分析、识别新的转录本和进行差异表达基因（differentially expressed genes，DEGs）鉴定等（Jeukenset al.，2010），还可以进行单核苷酸多态性的鉴定，开发大量简单重复序列标记和识别可变剪接（Garget al.，2011；Gaoet al.，2012）。Malenfant等（2013）通过构建白尾鹿Odocoileus virginianus血液转录组识别单核苷酸多态性，共鉴定出66 596 个SNP，丰富了其遗传资源信息；Jia 等（2016）对梅花鹿Cervus nippon进行了4个发育阶段的多组织转录组分析，揭示了其生命周期中最复杂、最重要的阶段是青少年过渡期；Yang 等（2016）对马鹿Cervus canadensis鹿茸组织进行从头转录组组装和分析，揭示了鹿茸的再生机制。为了解豚鹿基因表达的基本特征，本研究采用RNA-seq技术对成都动物园圈养的16 只豚鹿（6 雌10 雄）的血液进行转录组测序、从头组装，并进行基因功能注释、代谢通路分析以及雌雄差异表达基因分析，以期丰富豚鹿的基因注释信息并为其遗传多样性保护提供潜在的理论依据。

1 材料和方法

1.1 实验材料

16 份豚鹿血液样品（10 份雄性和6 份雌性）于2017—2018年采自成都动物园圈养种群，−80 ℃保存，送诺禾致源公司（北京，中国）使用Illumina HiSeqTM 2000 测序仪进行150 bp 双向测序，获得raw reads。

1.2 数据过滤及从头组装

测序获得的raw reads 进一步进行质量控制。当reads中未知碱基（N）的含量超过该条reads碱基数的10%以及reads 中的低质量碱基数超过该条reads 碱基数的50%时，去除此对reads，再使用Trimmomatic（Bolgeret al.，2014）去除含有接头污染的reads，最后得到高质量的clean reads。使用Trinity（Haaset al.，2013）将 clean reads 进行从头组装，Trinity 组装获得的转录组，一般含有相似的冗余序列，再使用Cd-hit-est（Li & Godzik，2006）去除冗余序列，得到的非冗余序列用于后续分析。

1.3 功能注释

为了获得转录本的注释信息，使用Blastx 将所有的转录本比对到 NR、GO、COG、Swiss-Prot 和KEGG数据库（E-value≤10−5）。为了解基因参与的生物学通路，所有的转录本上传到在线的KEGG 自动注释服务器进行注释（Moriyaet al.，2007）。

1.4 基因表达量分析及DEGs筛选

使用RSEM 1.2.2（Li& Dewey，2011）对每个样本基因表达量水平进行估计。表达量的结果使用FPKM 值表示。使用 DEseq2（Loveet al.，2014）进行不同样本之间的DEGs分析，筛选标准为log2|FC|>1且P<0.05。获得DEGs 集后，使用Bioconductor R语言包ClusterProfiler 3.8.1（Yuet al.，2012）对其进行GO 功能富集和KEGG 通路注释，以进一步评估DEGs的功能。

2 结果

2.1 转录组测序和组装

16 个血液样本提取RNA，利用Illumina HiSeq平台进行高通量测序，一共获得了118.6 G raw reads，包括 51 068 480 条 raw reads，去除被接头污染和低质量的序列后，得到了48 472 722 条共112.6 G clean reads。使用Trinity 将过滤后的clean reads进行从头组装，一共得到615 548条unigenes，长度为200~18 935 bp，平均长度为600.48 bp，Con‑tig N50 为803 bp，GC 含量平均值为46.89%。大多数转录本的长度为200~400 bp（260 743条），1 000~2 000 bp（50 248 条）的次之。由此可见，拼接组装后的转录组质量良好，可用于进一步的分析（表1）。

表1 豚鹿转录本组装结果统计Table 1 Summary of transcriptome assembly of Axis porcinus

2.2 转录本注释

与NR 的比对结果显示，有118 326 条转录本（37.93%）获得注释信息；在Swiss-Prot 中注释到的转录本为96 171 条（30.83%）；COG 中注释到的转录本为39 077条（12.53%）；在KEGG中获得注释信息的转录本共42 428 条（6.94%）。在四大数据库中均获得注释信息的转录本有6 307条（图1：A）。

GO分类显示15 971条转录本被分配到了生物过程（8 604 条）、分子功能（10 594 条）和细胞组分（5 855 条）。在生物过程中，细胞过程和代谢过程占主导；在分子功能中，转录本集中在连接和催化活性；在细胞组分中，转录本数量最丰富的类别是细胞、细胞部分和细胞器（图1：B）。

COG 中共有39 077条转录本被分配到25个功能类别中，数量最多的类别是一般性功能预测，其次是信号转导机制、氨基酸运输与代谢，而分布最少的是染色质的结构和动力学、细胞外结构和核结构（图1：C）。

KEGG 中注释的42 428 条转录本被划分到502 个代谢通路中，参与信号转导通路（ko09132）的最多，共有7 144 条；其次是病毒感染（ko09172，4 857 条）和与免疫系统相关的通路（ko09151，4 064 条）（图1：D）。根据免疫系统分类，免疫相关基因被映射到21 个免疫通路：涉及FcγR介导的吞噬作用（ko04666）的最多（1 894条），其次是自然杀伤细胞介导的细胞毒性（ko04650，1 880条）和B细胞受体信号传导途径（ko04662，1 823条）。

图1 转录本与数据库同源性比对注释结果Fig. 1 Results of homology search of transcripts against database

2.3 血液表达量最丰富的基因

根据NR数据库的比对并使用FPKM 值统计转录本表达水平，表达量最高的前10 个转录本主要包括主要组织相容性复合体（MHC）Ⅰ类抗原、几种核糖体蛋白和细胞色素氧化酶等基因，其中，HBA基因也是血液高表达的基因，编码了血红蛋白的α亚基（表2）。

表2 豚鹿血液转录组中表达量前10个转录本Table 2 The top 10 most abundant transcripts in the blood transcriptome of Axis porcinus

2.4 雌雄DEGs与富集分析

共鉴定到1 793 个性别DEGs。雄性显著上调的有1 781 个，显著下调的有14 个，显著下调的仅富集到1 条KEGG 信号通路：调节干细胞多能性的信号通路（ko04550），参与的基因为赖氨酸乙酰转移酶6A（KAT6A）。雄性显著上调的DEGs 的GO 富集显示，在分子功能类别中，DEGs主要富集在分子功能调节剂（GO：0098772）和酶调节活动（GO：0030234）；与细胞成分有关的DEGs 主要富集在与细胞骨架部分（GO：0044430）和微管组织中心（GO：0005815）有关的条目；在生物过程中，DEGs主要富集在免疫系统过程（GO：0002376）和囊泡介导转运（GO：0016192）（图2：A）。KEGG 富集显示，显著上调的DEGs 主要富集在血管内皮生长因子（VEGF）信号通路（ko04370）、趋化因子信号通路（ko04062）、溶酶体（ko04142）、mTOR 信号通路（ko04150）、破骨细胞分化（ko04380）、NOD 样受体信号通路（ko04621）和甘油磷脂代谢通路（ko00564）等通路上（图2：B）。其中，NOD 样受体信号通路、mTOR 信号通路、甘油磷脂代谢等通路含有部分免疫及性别相关基因，如参与免疫防御过程的TNFRSF25基因、与精子发育相关的精子鞭毛蛋白1（Spef1）基因，此外，与精子鞭毛形态形成有关的基因（Dnah17）以及睾丸表达基因（TEX264）也是雄性个体中上调的DEGs。这些雄性中高表达的DEGs可能与其生理特性相关。

图2 雄性豚鹿上调差异表达基因富集Fig. 2 Enrichment of up-regulated differentially expressed genes in male Axis porcinus

在雄性中显著高表达的DEGs 中，还发现了一些过去被报道参与鹿角生长发育相关的基因，包括与癌细胞生长有关的基因Rela、ADAMTS2，与神经嵴生长分化有关的基因OLIG1、Snai1，与软骨组织分化有关的基因PTH1R、TGFBI、MMP9，这些基因的表达量均显著高于雌性（图3）。

图3 部分差异表达基因在豚鹿中的FPKM 表达量Fig. 3 FPKM values of some differentially expressed genes of Axis porcinus

3 讨论

目前关于豚鹿的基因组虽已有报道（Wanget al.，2019），但该基因组缺乏基因注释信息。转录组的从头组装能够为基因表达提供重要信息，对基因组的补充注释具有重要意义。本研究基于16只豚鹿个体进行血液RNA-seq 测序，从头组装了一个血液转录组，N50 为803 bp，组装质量较好。基于不同的数据库对转录本进行注释，结果显示，仅有37.93%的序列获得了注释信息，低于水牛Bubalus bubalis（77.91%）（Denget al.，2016）、山 羊Capra hircus（45.41%）（Liuet al.，2013）等的注释率。这可能是由于参与组装的豚鹿个体多，测序的转录本中含有较多短序列；其次，注释信息不完善或可能存在豚鹿特有的新基因等。

豚鹿血液转录组大量表达信号转导、病毒感染和免疫系统相关基因，这与血液的重要生理功能密切相关。KEGG 注释显示许多转录本注释到与免疫相关的通路，这与大熊猫Ailuropoda melano⁃leuca（Duet al.，2015）、水牛（Denget al.，2016）等的结果一致，包括参与FcγR 介导的吞噬作用的大量基因。吞噬作用是重要的免疫机制之一，巨噬细胞在FcγR 受体作用下识别外源性感染物，触发特有免疫应答反应（Park，2003）。在表达量前10 的基因中，表达水平最高的是MHC Ⅰ类抗原复合体基因，MHC 分子在免疫应答过程中参与抗原识别，是动物免疫系统的重要组成部分（Kasahara，1999）。

通过比较雌雄性的血液转录组，雄性中显著上调的DEGs 有1 781 个，远多于显著下调的基因数量（14 个）。其中，雄性显著高表达的基因包括一些与精子形成、睾丸功能相关基因，如Dnah17与精子鞭毛发育有关（Shaet al.，2020）。此外，在KEGG 通路富集结果中，雄性上调的DEGs 富集到mTOR 信号通路，该通路是哺乳动物体内与细胞凋亡、细胞生长、免疫调节以及疾病发生等过程相关的重要通路（El Hianiet al.，2019），Mok等（2013）证实mTOR 信号通路对睾丸的发育调控具有重要作用，并且参与了精原细胞的分裂增殖。参与该通路的基因高表达可能与雄性个体生殖发育有关。雄性显著高表达的基因中还涉及许多免疫相关基因，如ZC3H12A基因的缺失与自身免疫性疾病有关（Cifuenteset al.，2010）。TNFRSF25参与免疫防御和炎症过程（Jaeckleet al.，2020），PSMB8编码免疫蛋白酶体的催化亚基，在人类白细胞抗原Ⅰ类呈递系统中起关键作用，与自身炎症性疾病相关（Kitamuraet al.，2011），肿瘤坏死因子（TNF）是免疫的中枢调节剂，通过介导细胞存活和细胞死亡诱导信号，参与免疫系统的发育和正常运作（Web‑ster & Vucic，2020）。而雌性只鉴定到1 个高表达的基因（KAT6A）与免疫相关，该基因的上调对于单核细胞向巨噬细胞分化和促炎细胞因子的产生很重要（Wiesel-Motiuk & Assaraf，2020）。这些结果表明，豚鹿在自身免疫方面存在显著性别差异，这可由哺乳动物的免疫系统具有广泛的性别二态性来解释（Gal-Ozet al.，2019），不同的因素包括遗传因素和激素介质都可以解释免疫反应中基于性别的差异（Oertelt-Prigione，2012）。

雌雄豚鹿的DEGs 还揭示了与鹿角生长发育相关的基因。KEGG 富集发现，雄鹿中上调的DEGs显著富集到血管内皮生长因子信号通路和破骨细胞分化等通路上。破骨细胞分化的调节对正常的骨重塑至关重要（Sharmaet al.，2006）。血管内皮生长因子（VEGF）具有促进血管形成和骨形成的作用，刺激内皮细胞的分裂增殖，能促进新生血管的形成并增加血管通透性（Apteet al.，2019）。鹿角的生长需要血管大量供血，因此富集到VEGF信号通路的DEGs 可能与鹿角的生长形成有关。通过基因的数据库比对信息，在雄性上调的DEGs中，鉴定了可能与鹿角生长相关的基因：Rela、ADAMTS2、OLIG1、Snai1、PTH1R、TGFBI、MMP9。Rela是一种原癌基因（Lu & Yarbrough，2015），Wang等（2019）的研究认为鹿角生长和肿瘤细胞发生具有相似的发育模式，鹿角组织的快速生长特性与肿瘤发生途径有关。因此该原癌基因在雄性中的高表达可能与鹿角的快速生长有关，以满足鹿角快速生长的代谢需求。ADAMTS2基因是ADAMTS家族成员，Wang 等（2019）指出ADAMTS家族成员中的ADAMTS2、ADAMTS4、ADAMTS12、ADAMTS14、ADAMTS17和ADAMTS18是鹿角的特异性基因，在鹿角中的表达量很高。它们通过控制细胞外基质的结构和功能以及调节肿瘤微环境来抑制癌细胞的生长（Jinet al.，2007），该基因在雄鹿中的高表达可能是为了防止与鹿角生长相关的细胞因快速增殖而发生癌变。OLIG1基因与神经嵴分化途径有关（Betancuret al.，2010），神经嵴细胞具有迅速增殖和分化为软骨和神经细胞的潜力，Wang 等（2019）发现鹿角基因都是从神经嵴干细胞发育而来，且Carlson 等（2016）的研究表明，OLIG1基因突变导致了牛的无角性。因此猜测OLIG1基因可能也在豚鹿鹿角的进化和发育过程中起着至关重要的作用。Snai1基因是与神经嵴细胞迁移相关的基因，被证实在胎儿羊角中高表达（Wanget al.，2019），PTH1R和TGFBI基因在鹿角软骨组织中高表达，并且能调节软骨细胞和破骨细胞的增殖和分化（Faucheuxet al.，2002，2004），MMP9基因在鹿角软骨细胞中高度表达，参与基质降解（Faucheuxet al.，2001）。这些与鹿角发育相关DEGs 的表达量均在雄鹿中上调，表明其可能在雄性豚鹿鹿角形成和发育中发挥着重要作用，但它们在鹿角生长过程中的具体调控机理还有待进一步挖掘探讨与功能验证。本研究结果不仅为丰富豚鹿转录表达提供了基础数据，也为揭示偶蹄目和鹿科动物性别差异和鹿角形成机制奠定了基础。