基于RNA-seq 对复合益生菌作用于肉鸡回肠中新转录本预测及可变剪接、SNP 分析

曹威荣， 王燕飞， 张若男， 贾 浩， 刘 璇， 张利环

（1.朔州职业技术学院，山西 朔州 036002；2.山西农业大学生命科学学院，山西 晋中 030801）

黄羽肉鸡抗病力和对劣质饲料的适应性强，易饲养，具有较高的繁殖性能，肉质风味与口感也较佳。我国黄羽肉鸡种类繁多，其中黄麻肉鸡体型圆润、整体发育匀称、鸡皮较薄、皮下脂肪适中，肉中富含较高蛋白和维生素A，具早熟易肥、肉质紧实滑嫩等特点，因此在全国被广泛养殖（阿卜杜热合曼·达尼亚尔，2018）。 食用动物摄取的抗生素，以及产生的食物中的抗生素残留， 已被公认为人类耐药性迅速增加的主要原因之一 （Davies 等，2010）。 滥用抗生素饲喂食用动物，会直接选择对抗生素有抗药性的微生物， 并把食用动物的系统转变成抗生素抗药性基因的贮存库（Claesson 等，2012）。因此，减少食用抗生素的摄入量，消除食用动物农业中的抗生素残留， 已成为食品安全和公共卫生的首要任务。

目前，在养鸡业中，饲料中抗生素的替代品包括纤维降解酶、益生菌、益生元、共生元和植物生物元素等（Sopkova 等，2017）。 其中益生菌生产成本低，在不同种类的宿主动物中应用范围广，具有优势（Gaggia 等，2010）。益生菌被定义为在机体内定植一段时间后可改善肠道微生态并对宿主动物产生积极影响的活性微生物 （Sanders 等，2008）。干酪乳杆菌、 嗜酸乳杆菌和双歧杆菌一起被称为“健康三益菌”。 乳酸菌是一种非致病性和毒素性的革兰氏阳性发酵菌，无芽孢与鞭毛，菌落白色粗糙，区别于其他厌氧菌，大多数可以在氧气存在下生长。双歧杆菌是定居于人胃肠道的微生物之一，虽不属于乳酸菌， 但可以产生乳酸。 在饮用含有1×109CFU 的干酪乳杆菌溶液时， 益生菌因其特性对胃液、 水解酶和胆汁酸具有高耐受性并可在胃肠道中存活并停留3 d（Radicioni 等，2019）。 乳酸菌在肠道中因可分泌诸如细菌素的抗菌物质而具有抗炎性和抑菌特性，产生的酸可降低肠道pH并防止产气荚膜梭菌的增殖，有利于肠道健康。

RNA-seq 数据的高分辨率不仅可以探索基因表达， 还可以预测动物中的可变剪接事件（Li等，2019）。 RNA-seq 可以反映出mRNA、smallRNA、noncodingRNA 等或者其中一些的表达水平，有助于查看突变/SNP 和基因表达随时间的变化（Maher 等，2009），还可以查看miRNA、tRNA 和核糖体分析等（Ingolia 等，2012）。SNP 是一种基于核苷酸变异性的有效遗传标记， 在动物遗传学研究和育种中有着广泛的应用（Li 等，2019）。 可变剪接是真核生物中的一种正常现象， 其带来了重要的蛋白质多样性并允许基因产生不同的mRNA转录本， 这些转录本被翻译成不同的蛋白质（Li等，2020）。

本试验对肉鸡回肠进行了新转录本的发掘和注释、可变剪接和SNP 分析，以期为益生菌饲喂肉鸡新基因的挖掘提供数据基础， 进而更好的探究其分子机制。

1 材料与方法

1.1 试验动物及样品采集 干酪乳杆菌、嗜酸乳杆菌与双歧杆菌冻干菌粉有效活菌数为1×1010CFU/g， 冻干菌粉购自陕西瑞茂生物科技有限责任公司，复合菌制剂（干酪乳杆菌：嗜酸乳杆菌：双歧杆菌=1:1:2）（Zhang 等，2022）。200 只1 日龄黄麻肉鸡（公母各半，购自太谷县田建胜养殖场）随机分为2 组，每组5 个重复，每个重复20 只。 对照组正常饮水，益生菌组饮用水中补充1%复合菌制剂，试验分为生长前期（1 ～21 d）、后期（22 ～42 d）两个阶段。 基础日粮配制参照NRC（1994）和《NY/T 33-2004 鸡饲养标准》（表1）。 42 日龄时进行肉鸡屠宰试验，取回肠后用生理盐水将其冲洗干净，分装后液氮速冻，并于-80 ℃保存。

表1 基础日粮组成和营养水平

1.2 RNA-Seq 测序 从对照组和益生菌组随机抽取3 个样本 （编号为C1、C2、C3 和P1、P2、P3）送至上海美吉生物医药科技有限公司，使用Illumina Novaseq 6000 测序平台进行转录组测序（RNA sequencing，RNA-Seq）。

1.3 新转录本的发掘 使用Cufflinks 软件对Mapped Reads 进行拼接，并与原有的基因组注释信息进行比较，寻找原来未被注释的转录区，发掘该物种的新转录本和新基因。

1.4 新转录本功能注释 使用BLAST 软件将发掘的新转录本分别与KEGG（http://www.genome.jp/kegg/）、COGs （http://www.ncbi.nlm.nih.gov/COG/）、GO （http://www.geneontology.org/）、NR（ftp://ftp.ncbi.nlm.nih.gov/blast/db/）、Swiss -Prot （ftp://ftp.uniprot.org/pub/databases/uniprot/current_release）、Pfam（http://pfam.xfam.org/）等数据库进行比对来获得测序所得肉鸡新转录本的注释信息。

1.5 可变剪接分析及差异可变剪接基因分析使用rMATS（http://rnaseq-mats.sourceforge.net/index.html）对RNA-seq 数据进行可变剪接分析。利用GO、KEGG 数据库进行差异可变剪接基因富集分析。

1.6 SNP 分析 使用Sentieon 进行SNP 分析，得到SNP 信息。

2 结果与分析

2.1 测序数据质控结果及与参考基因组比对分析 如表2 所示，样本过滤后的Clean reads 数占Raw reads 的比例为98.26% ～98.91%， 且Clean bases 占Raw bases 的比例为96.16% ～97.48%；测序碱基平均错误率均为0.026%以下；碱基识别准确率在99%以上的碱基所占百分比均在97.51%以上， 碱基识别准确率在99.9%以上的碱基所占百分比均在93.27%以上，GC 含量在样本中所占比例为49.09% ～51.12%，没有分离现象，说明测序产出数据质量较高，符合测序要求。

表2 序列质量和比对信息汇总统计

测序数据与参考基因组的比对率为89.04% ～91.48%， 说明不存在污染且比对基因组正确；在参考序列上有唯一位置的比对率为86.88% ～89.35%， 在参考序列上有多个位置的比对率为2.13% ～2.39%，说明基因比对效果较好，且基因覆盖均一度较好。

2.2 新转录本的发掘及注释 如图1 所示，将测序数据与参考基因组注释信息进行比对， 共发现1047 个新基因， 注释到GO、KEGG、COG、NR、Swiss-Prot、Pfam 数据库中的基因分别有276、196、323、567、250、247 个。

图1 各数据库注释的新基因数

如表3 所示， 将高通量测序数据与参考基因组注释信息进行比较，共发现20176 个新转录本，注释到GO、KEGG、COG、NR、Swiss-Prot、Pfam 数据库的新转录本数分别为13739、13488、17693、18764、17499、16655。

表3 各数据库注释的新转录本数

如图2 所示， 转录组测序得到的新基因共276 个基因在GO 数据库中得到归类注释。 分子功能主要集中在结合性、催化活性等多个分类中；在细胞组分分类中，主要集中于生物膜组分、细胞成分、细胞器等中；在生物学过程分类中，主要集中在细胞过程、代谢过程、生物调控等中。

图2 新基因GO 分类注释图

2.3 肉鸡可变剪接分析及差异可变剪接基因分析 如图3 所示， 其中SE 占可变剪接类型的比例最高， 达到70.11% ～70.28%，A5SS、A3SS、MXE、RI 分别占各个阶段可变剪接事件总量的6.80% ～7.75%、1.04% ～1.19%、5.80% ～6.12%、4.49% ～5.18%。

图3 可变剪接事件统计图

如图4 所示，共筛选到1131 个显著的差异剪接基因， 在SE、A5SS、A3SS、MXE、RI 事件中鉴定到的差异剪接基因数分别为588、139、176、136、92。

图4 差异剪切基因统计图

如图5 所示，在SE 中外显子排斥有327 个，外显子包含有410 个； 在A3SS 中外显子排斥有99 个，外显子包含有101 个；在A5SS 中外显子排斥有72 个，外显子包含有87 个；在MXE 中外显子排斥有138 个，外显子包含有117 个；在RI 中外显子排斥有55 个，外显子包含有58 个。

图5 可变剪接事件数目外显子纳入和排斥统计图

如图6 所示， 对获得的差异剪接基因进行GO 富集分析，图中展示了前20 个GO 条目，富集到了与代谢和免疫相关的GO 条目， 如大分子代谢过程的调控、 磷脂酰肌醇磷酸4-磷酸酶活性、细胞代谢过程、免疫反应的正调控、大分子代谢过程的正调控。

图6 差异剪切基因GO 富集柱状图

如图7 所示，差异可变剪接基因KEGG 富集到了以下通路， 胞吞作用、MAPK 信号通路、Rap1信号通路、GnRH 信号通路、 黏附连接、TNF 信号通路、磷脂酰肌醇信号系统、ErbB 信号通路等。

图7 差异剪切基因KEGG 富集气泡图

2.4 SNP 分析 由表4 可知， 转换SNP 所占百分比为73.53% ～73.94%； 颠换型SNP 所占百分比为26.05% ～26.47%， 在各样品中大多数碱基取代均为转换大于颠换。

表4 SNP 类型统计

3 讨论

转录组学研究的是在特定发育阶段或生理条件下细胞RNA 水平上基因的表达和变化情况，进而查找探究基因的功能， 它在新转录本预测方面有广泛应用（Chen 等，2011），还用于小菜蛾、帝王蝶等昆虫新测序基因组的基因注释 （You 等，2013；Zhan 等，2011）。 Izadnia 等 （2019） 利用RNA-seq 技术鉴定罗斯708 和伊斯法罕鸡两种鸡肝脏与剩余采食量（RFI）相关的关键基因和新的转录体。Liu 等（2020）在鹧鸪鸡中37.6%的蛋白质编码基因中发现了新的转录体。 Tang 等（2020）对来自Valgus-varus 畸形哈伯德肉鸡和哈伯德肉鸡的6 个鸡脾脏样本进行了转录组测序发现了7943 个新的转录本。 本研究共发现20176 个新转录本，1047 个新基因， 新基因中共276 个基因在GO 数据库中得到归类注释。

目前， 转录组测序已成为分子生物学研究的常用工具， 广泛应用于检测新的转录本（Huang等，2010）、差异基因筛选（Filichkin 等，2010）、可变剪接 （Lan 等，2014） 以及SNP 筛查（Martinez等，2017）等相关研究。 真核生物中存在着大量的可变剪接（Li 等，2016），动物细胞中也存在着组织特异性（Montelli 等，2015）。Hu 等（2021）对北京鸭胸肌和腿肌进行转录组测序分析发现胸肌肌肉组织的可变剪接数量为32256 ～45707，腿部肌肉组织的可变剪接数量在不同时间点为35841 ～47889。Li 等（2015）发现，A3SS、A5SS、RI、SE 是20周龄青年鸡和30 周龄产蛋鸡肝脏中主要剪切的模式， 四类剪切事件占总数的96.5%，MEX 仅占3.5%。 Hu 等（2021）预测了汉中麻鸭不同生长阶段的骨骼肌可变剪接事件，在17、21、27 d 和6 个月大时的骨骼肌中可变剪接最多 （＞10000）。Daines 等（2010）利用RNA-Seq 数据与黑腹果蝇基因组比对，识别了至少4618 个基因（31%）发生了可变剪接， 确定了7776 个潜在的新剪切事件，最常见的带注释的事件是外显子跳跃（SE）。 本研究中SE 的可变剪接类型的比例最高， 共筛选到1131 个显著的差异剪接基因，并对差异可变剪接基因进行GO 富集分析， 主要富集到代谢和免疫GO term。

SNP 是指在基因组上由单个核苷酸的变异所引起的DNA 序列多态性，其数量很多，多态性丰富， 主要发生的有2 种， 即转换和颠换。 黄育雯（2020） 在安格斯牛下丘脑和垂体研究中发现，碱基转换的可能性超过了颠换。 Zhang 等（2019）检测了大横肉鸡MSTN 基因的外显子序列， 共检测到5 个SNPs。 张蕾等（2019）利用转录组测序技术对6 个高邮鸭卵巢组织进行研究发现， 转换类型总数极显著高于颠换类型总数。 姚娜等（2012）发现， 西藏班戈绒山羊和辽宁绒山羊山羊两品种间有5391 个SNP 位点。本研究在6 个样品上分别得到了可能的SNP 位点，SNP 转换类型总数大于颠换类型总数。

4 结论

本试验对肉鸡的回肠进行了RNA-seq 测序分析，与参考基因组和参考基因比对发现了1047个新基因，20176 个新转录本； 共筛选到1131 个显著的差异剪接基因， 在SE、A5SS、A3SS、MXE、RI 事件中鉴定到的差异剪接基因数分别为588、139、176、136、92；SNP 中转换类型总数大于颠换类型总数。 本研究为发掘复合益生菌作用于肉鸡回肠中的可变剪接事件提供基础， 为进一步了解复合益生菌作用于肉鸡回肠中新转录本的发现和SNP 位点提供了数据支撑。