易 立，仝明薇，程悦宁，程世鹏

(中国农业科学院特产研究所 特种经济动物分子生物学国家重点实验室，吉林 长春 130112)

乌苏里貉是貉属的一种珍贵毛皮动物，现阶段我国广泛饲养的乌苏里貉亚种都是地产品种[1]。目前，我国貉的免疫相关分子信息挖掘亟待开展，如Toll 样受体、干扰素、细胞因子及受体等貉免疫分子信息大部分尚未解析，采用传统的依据同源基因(相似动物)设计引物，克隆测序显得较为繁琐。采用2 代测序通过一次测序就能获得大量序列信息，结合生物信息学比对，获得全基因组序列信息及注释信息[2-3]。

转录组是一个细胞中的所有转录本信息，包括转录本的数量、表达动态、转录本序列信息。由于转录组(mRNA 测序)几乎所有序列均为外显子，去除了内含子等无效信息，所以更为精简有用。乌苏里貉尚无全基因组序列，因此，建立貉外周血单核细胞(PBMCs)转录组测序及分析，对于研究其免疫系统特征及免疫机制具有重要意义。

本研究利用Illumina 平台首次对貉PBMCs 转录组进行序列分析，并对RNA-seq 数据进行拼接和注释，为后续开展貉进化研究及免疫机制研究提供参考，为貉全基因组测序项目提供较好的注释基础。

1 材料和方法

1.1 主要实验材料 外周血样品采集自10 只配种后公貉，每只约20 mL，共计200 mL；红细胞裂解液购自康为世纪公司；TRIzol Reagent 购自Invitrogen公司；PolyATtract®mRNA Isolation System III with Magnetic Stand 试剂盒购自Promega 公司；TransScript®One-Step gDNA Removal and cDNA Synthesis SuperMix试剂盒购自TransGen Biotech 公司；QIAquick PCR Purification 试剂盒购自QIAGEN 公司。

1.2 貉PBMCs 总RNA 提取 按照红细胞裂解液说明书方法，处理貉外周血并收获PBMCs，参照TRIzol Reagent 说明书提取总RNA。利用1.2 %琼脂糖凝胶电泳及分光光度计检测确定总RNA 完整性和纯度质量。

1.3 cDNA 文库的构建和Illumina 测序 提取样品总RNA，采用带有Oligo(dT)的磁珠富集mRNA；在制备的mRNA 中加入fragmentation buffer 使其片段化为短片段，以其为模板，利用六碱基随机引物合成cDNA 第一链，并加入缓冲液、dNTPs、RNase H和DNA polymerase I 合成cDNA 第二链，经过QIAquick PCR Purification 试剂盒纯化并加EB 缓冲液洗脱经末端修复、加碱基A，加测序接头，回收目的片段后进行PCR 扩增，从而完成整个文库制备工作。构建的文库利用Illumina HiSeqTM2500 进行测序，测序策略为PE100(Paired-end 技术，双向测序100 bp)测序[4]。

1.4 质量控制和序列的拼接 测序得到的某些原始序列含有带接头的、低质量的序列，为保证信息分析质量，对原始序列进行过滤，得到clean reads，后续分析均基于clean reads。数据处理步骤包括：1)去除含N 比例大于10 %的read；2)去除超过50 %的碱基质量值低于5 的read；3)去除接头污染。对质量控制前后整体序列的质量进行统计，质量控制前后reads 的数量和长度通过perl 脚本进行统计。

序列的de novo 拼接采用基于构建De Bruijn 图方法的Trintiy[5]，Trinity 能够高效的构建de novo 转录组。Trinity 利用de Bruijn 图论的原理，针对转录本具有可变剪接的特点，组装出全长转录本。

1.5 转录物功能注释及分类 利用Trinotate 对开放阅读框(ORF)和contigs 进行功能注释，利用Uniprot database、RNAMMER、eggNOG、GO、KEGG 对预测的序列进行注释。由于目前无貉基因组信息，因此本研究中的序列拼接参考BLASTX、COG(Clusters of orthologous groups of proteins)、GO(Gene ontology)数据库的比对结果，利用BLASTX 将组装序列与COG 直系同源簇数据库进行比对分析，获得貉表达基因的COG 功能注释及其COG 功能分类。最后分别使用Blast2GO 软件和WEGO 工具对转录物序列进行GO 功能注释及功能分类统计。

1.6 ORF 预测 对组装的转录本利用TransDecoder鉴定编码区域，并对长度大于300 bp 的序列进行best ORF 查找。

1.7 基因表达分析 使用bowtie 软件，并按照每个样品各自特异的测序接头信息把对应的reads mapping 到已拼接的转录物上，以转录物序列中的reads 数表示基因的表达丰度；分析转录组整体表达情况，利用组装的转录本，采用RSEM(RNA-Seq by Expectation-Maximization)软件估计每个样品表达量。

2 结果

2.1 总RNA 质量检测 按常规方法提取貉PBMCs样品总RNA 并对其浓度进行检测，其OD260nm/OD280nm=2.05，浓度为212 ng/μL，表明提取的RNA质量较好，可以满足正常的实验要求。

2.2 测序数据初步分析 Illumina 测序共得到3.1 GB 的数据，共32 245 804 条未加工数据(Raw reads)，经过质量预处理后共得到有效数据(Valid reads)28 797 350 条。经过Trinity 拼接去重复后，共得到118 868 条长度大于200 bp 的contigs，平均长度为525.53 nt，N50 值为674，最长的contig 为8 982 bp，最短的为305 bp，contigs 的长度统计分析见表1。

2.3 转录物功能注释及分类 为从整体上了解转录物序列信息，首先利用Uniprot database 数据库对上述转录物数据进行BLASTX 比对。随后利用COGs、GO 及功能注释分类体系对具有同源对比信息的转录物进行比对和功能注释(E-value<1e-05)。在COG功能分类体系中，共获得8 261 个COG 功能注释，涉及24 个COG 功能类别(图1)。其中，一般功能基因(General function prediction only)的转录物比例最大；其次为蛋白质翻译后修饰与转运及分子伴侣相关基因(Posttranslational modification，protein turnover and chaperones)。在该转录组中，涉及貉生长发育相关的功能定义主要包括：氨基酸转运与代谢(E)、碳水化合物转运与代谢(G)、核酸转运与代谢(E)和脂质转运与代谢(I)等物质代谢过程、翻译，核糖体结构与生物合成(J)、蛋白质翻译后修饰与转运，分子伴侣(O)、信号转导机制(T)、无机离子转运与代谢(P)、辅酶转运与代谢(G)、次级代谢产物生物合成、运输与分解代谢(Q)及防御机理(V)等多个生理生化过程。其中信号转导机制(T)及防御机理(V)功能定义可能与免疫信号相关。S 为功能未知基因，从中可以发掘新基因。

图1 貉PBMCs 转录组24 种COG 功能分类Fig.1 Twenty-four types of COG function classification of the PBMCs transcriptome of raccoon dog

对所有转录本进行GO 功能分类(图2)。转录物的GO 注释归为生物学过程、细胞组分和分子功能3 大功能，分别包含了23、14 和17 个功能亚类。在生物学过程功能类型中，主要的生命过程分别是细胞过程(Cellular process)和代谢过程(Metabolic process)，细胞杀伤(Cell killing)及节律过程(Rhythmic process)比例最低。在细胞组分功能类型中，细胞(Cell)和细胞部分(Cell part)所含比例最高；而胞外区部分(Extracellular region part)、合胞体(Symplast)、神经元(Synapase)、神经元组成(Synapase part)、病毒(Virion)、病毒组成(Virion part)均未检测到。在分子功能类型中，蛋白结合(Binding)和转录调节活性(Transcription regulator)所含比例最高；抗氧化活性(Antioxidant)、辅助运输蛋白活性(Auxiliary transport protein)、蛋白标签(Protein tag)、调节酶活性(Proteasome regulator)等均未检测到。

图2 貉PBMCs 转录组的54 种GO 功能分类Fig.2 54 types of GO classification of the PBMCs transcriptome of raccoon dog

通过KEGG 数据库，本转录组共有16 650 条能够在KEGG 数据库中得到注释，并分别定位到259条相应的通路上。其中，T-细胞受体信号通路、B-细胞受体信号通路、补体-血凝联级反应通路、吞噬体通路、自然杀伤细胞介导的细胞毒性通路、产生IgA 的肠道免疫网络通路、趋化因子信号通路、TLR 受体信号通路、NLR 受体信号通路、RIR 受体信号通路、溶酶体通路等通路与貉免疫与抗病相关。其中Toll 样受体信号途径(Toll-Like receptor signaling pathway)见图3，在该信号途径中，绝大部分的分子均能够表达。

2.4 ORF 预测及基因表达分析 为获得具有完整ORF 的重要功能基因序列，本研究选用TransDecoder对转录物进行分析鉴定编码区域，共获得118 868条ORF 信息。依据序列中包含的reads 个数(表达量)，对转录物表达丰度进行分析。转录物reads 个数分布范围较广，FPKM+1 取log10 为表达量，制作表达量密度分布图(图4)，表达量最高的序列在logFPKM 值为0.6～0.8，表明此次测序的重复率较低。

图3 貉PBMCs 转录组Toll 样受体信号途径分子表达Fig.3 Toll-Like receptor signaling pathway of the PBMCs transcriptome of raccoon dog

图4 貉PBMCs 转录组表达量密度图Fig.4 Expression density of the PBMCs transcriptome of raccoon dog

2.5 本地BLAST 数据库的建立及初步应用 下载NCBI 本地BLAST 软件，将拼接的118 868 条ORF信息通过Makeblastdb 命令进行数据的初始化，建立貉PBMCs 转录组数据库。利用本地BLAST 命令，对其进行检测，以貉Toll-1 基因为例，进行本地BLAST，共搜寻到7 条匹对的组装后序列，相似率均为100 %。表明本研究建立的数据库可信。

3 讨论

本研究在文库构建中，采用RNA 样本片段化后再反转处理的方法，充分利用了RNA 对二价阳离子的敏感性及稳定性好、减少RNA 二级结构等优点，获得了更均匀的覆盖率和更全面的转录物信息[6]。采用双端测序(Paired-end)方法进行高通量测序，不但增加了测序的深度，而且提高了de novo 拼接的效率和准确性。

本研究利用Uniprot database 数据库对转录物数据进行BLASTX 比对，结果显示，有少数序列与其它物种蛋白序列并无匹配，分析原因如下：1)序列片段太短，无法获得同源性比对结果。随着序列长度的增加，能比对上的序列比例会明显增加。芝麻[7]、水稻[8]等物种的转录组序列分析中均出现过此现象。2)注释信息不完全。目前貉基因组学及转录组学研究较为空白，其生物信息数据库仍在不断的更新完善中，基因功能注释信息不全面造成了部分序列暂时无法获得功能注释信息。3)貉物种或其组织转录组的特异性问题。貉基因组学及功能基因组学研究相对落后，目前常选用犬基因组作为参考序列，但由于犬与貉在遗传进化上距离相对较远，在貉转录本中可能存在某些属特异性新基因，导致其同源序列难以被发现。

对转录本进行功能注释(E-value<1e-05)。在COG 及KEGG 功能分类体系中，有多种注释涉及貉免疫功能，因此，后续的研究重点将放在进一步分析每种注释中匹配的基因，发掘参与貉免疫应答和防治疾病过程中的关键基因。

根据研究物种是否有参考基因组信息，RNASeq 序列拼接的策略也分为基于参考基因组的序列拼接和de novo 序列拼接[9-10]，或当基因组信息不完整时，将二者结合起来进行序列拼接。De novo 策略使用的软件有Trinity 和Oases 等，其中Trinity 可以发现更多的全长转录本，其灵敏度甚至能够近乎于基于基因组信息的拼接。目前已经发表的利用de novo 拼接方法进行的非模式生物转录组研究涵盖了昆虫、植物等许多物种[11-12]。

本研究首次对貉PBMCs 进行了较全面的转录组测序，PBMCs 作为重要的免疫细胞，其转录本中可能包含更丰富的免疫相关基因信息，为后续开展貉进化研究、免疫机制研究提供参考并为貉全基因组测序项目提供较好的注释基础。

[1]朴厚坤.特种经济动物养殖业的现状与展望[J].吉林农业大学学报，1995，(1)：100-105.

[2]杜玲，刘刚，陆健，等.高通量测序技术的发展及其在生命科学中的应用[J].中国畜牧兽医，2014，12：109-116.

[3]岳桂东，高强，罗龙海，等.高通量测序技术在动植物研究领域中的应用[J].中国科学：生命科学，2012，02：107-124.

[4]Marioni J C,Mason C E,Mane S M,et al.RNA-seq:An assessment of technical reproducibility and comparison with gene expression arrays[J].Genome Res,2008,18:1509-1517.

[5]Haas B J,Papanicolaou A,Yassour M,et al.De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis[J].Nature Protocols,2013,8(8):1494-1512.

[6]Wang Z,Gerstein M,Snyder M.RNA-Seq:A revolutionary tool for transcriptomes[J].Nat Rev Genet,2009,10(1):57-63.

[7]Wen Liang-wei,Xiao Qiong-qi,Lin Hai-wang,et al.Characterization of the sesame(Sesamum indicum L.)global transcriptome using Illumina paired-end sequencing and development of ESTSSR markers[J].BMC Genomics,2011,12:451.

[8]Ting Ting-lu,Guo Jun-lu,Dan Lin-fan,et al.Function annotation of the rice transcriptome at single-nucleotide resolution by RNA-seq[J].Genome Res,2010,20:1238-1249.

[9]Birzele F,Schaub J,Rust W,et al.Into the unknown:Expression profiling without genome sequence information in CHO by next generation sequencing[J].Nucleic Acids Res,2010,38(12):3999-4010.

[10]Sultan M,Schulz M H,Richard H,et al.A global view of gene activity and alternative splicing by deep sequencing of the human transcriptome[J].Science,2008,5891:956-960.

[11]Bin Chen,Yu Juan-zhang,Zheng Bo-he,et al.De novo transcriptome sequencing and sequence analysis of the malaria vector Anopheles sinensis(Diptera:Culicidae)[J].Parasit Vectors,2014,7(1):314.

[12]Rastogi S,Meena S,Bhattacharya A,et al.De novo sequencing and comparative analysis of holy and sweet basil transcriptomes[J].BMC Genomics,2014,15:588.