摘要
为探究可变剪接影响下更为复杂的脂肪沉积调控网络，明晰可变剪接事件的发生对脂肪沉积调控机制的影响，本研究基于PacBio测序平台的第三代测序技术，对夷陵黄牛腹腔脂肪、皮下脂肪、肌间脂肪进行全长转录组测序分析。共发现15 445个基因检测到50 520个可变剪接，占牛全部基因的33.4%。对这些基因进行富集分析发现，83个GO条目显著富集，这些显著富集的通路可能参与脂肪沉积网络的调控，其中16个与脂质合成及代谢相关，27条KEGG通路显著富集，其中15条与脂质合成及代谢相关。使用KOG、KEGG、NR、Swiss Prot、GO数据库对测序序列进行注释，80 756 条序列共注释到 69 259 个基因，94 458条CDS序列及对应的pep序列。其中GO分析中15 039条序列富集到507个生物学过程条目，7 907条序列富集到214个细胞组分条目，30 672条序列富集到559个分子功能条目。KEGG数据分析富集到34条通路共49 710条序列，其中7 002条序列参与信号转导途径。以上结果表明，牛脂肪组织存在大量可变剪接事件，并且这些发生可变剪接事件的基因对脂肪沉积的调控发挥重要作用，这可从基因可变剪接角度为解析脂肪沉积的调控网络提供必要的资源和理论依据。
肌内脂肪的含量以及沉积分布影响牛肉的嫩度、系水力、剪切力值、风味和多汁性，适当提高牛肉肌内脂肪含量可改变牛肉品
可变剪接是剪接体对前体mRNA选择性的修饰，使同一基因可以表达不同成熟mRNA剪接异构体的现
基于PacBio 平台的全长转录组测序是第三代测序技术的代表，能直接测序得到RNA转录本的全长片段，可以规避短片段带来的缺失信息的影响，提供更完整的转录组信息，更好地进行差异表达基因分析及功能注
供试牛为湖北省宜昌地区、年龄在3岁、发育正常、体态良好的3头健康夷陵黄牛。屠宰后，使用消毒手术刀，分别在每头牛相同部位采集5 c
使用TRIZOL法提取夷陵黄牛脂肪组织总RNA。提取后的RNA需使用1%凝胶电泳检测样品纯净度，并经过Nano Drop仪器检测，确定总RNA的OD260/OD280值在 1.9~2.1、OD260/OD230值在1.4~1.8。并使用Qubit对RNA进行定量，Agilent 2100、Bioanalyzer 等进一步进行质量检测，提取的总 RNA质量满足要求后，低温保存，待后续进行文库构建。
构建文库所使用的RNA需要使用 Qiagen 试剂盒进行回收纯化预处理，以获得更高质量的 RNA。随后使用SMARTer
对建库测序得到的原始数据，首先使用SMRTLINK 5.1 软件对其进行质控，然后使用SNR（Signal Noise Ratio）筛选剔除低质量数据及接头序列，得到环形一致序列（circular consensus sequence， CCS）。去除含 5' primer、3' primer、poly A 片段，获得的全长转录组序列经过校正和去冗余后，可用于后续分析。
试验使用ICE （iterative isoform-clustering）中的algorithm 模块进行聚类，使用 DAGCon 得到一致性序列，并通过arrow 进行校正。将二代数据校正后的全长序列与参考基因组比对，去除融合基因及非注释数据后，将比对到参考基因组同一基因上的不同isoform进行提取，即为可变剪接。根据可变剪接不同的形成方式进行分类。
对采集的夷陵黄牛的皮下脂肪、腹腔脂肪、肌间脂肪样品使用PacBio平台第三代测序技术进行全长转录组测序，共获得19.81 G的原始数据，筛除接头和长度小于50 bp的原始离线数据后，得到13 014 954条Subreads，平均长度为 1 474 bp（Full Passes≥1 、最小预测准确度0.8），核验下机数据中的CCS（circular consensus sequence）的质量，共获得779 345个CCS reads，CCS总碱基数为1 812 898 181 bp，CCS Read 平均长度2 320.43 bp，平均 Full Pass 数量为11。通过验证5′引物、3′引物和poly-A尾的保留，获得高质量的FLNC（full-length nonchimeric reads）640 205条，占比 82.15%，平均长度为2 091 bp，以上结果说明全长转录组测序数据质量较好，可以用于后续试验。
利用来源于相同样品的二代测序数据对PacBio平台的第三代测序reads进行校
校正后的数据与参考基因组进行比较，得到与2 095个融合基因相关的2 768条reads。将剩余数据进一步过滤去冗余。使用Match Annot软件将比对后结果与注释信息进行比较，并将三代注释结果和原基因组结果进行合并，得到基因注释的 PB isoform 63 400 条，非基因区序列有36 480 条。将比对到参考基因组cDNA区域的序列进行检测，共识别到可变剪接事件50 520个。其中2 406个基因发生外显子跳跃、1 184个基因发生外显子替换、1 262个基因发生外显子接受、4 085个基因发生内含子保留、2 365个基因发生外显子位置替换、4 143个基因发生其他可变剪接事件（
| 类型Type | 剪接数目 Number of AS | 基因数 Number of gene |
|---|---|---|
| 外显子跳跃 Exon S | 4 677 | 2 406 |
| 外显子替换 Alt D | 1 939 | 1 184 |
| 外显子接受 Alt A | 2 138 | 1 262 |
| 内含子保留 Intron R | 14 606 | 4 085 |
| 外显子位置替换 Alt P | 6 490 | 2 365 |
| 其他 Others | 20 670 | 4 143 |
注： AS 统计根据测序注释结果与参考基因转录本比较得到。Note： AS statistics are based on the comparison of sequencing annotation results with reference gene transcripts.
将检测到可变剪接事件的基因进行GO富集分析，显示与83个条目显著相关（P<0.05），其中富集到分子功能类的GO条目有22个，细胞组分类的GO条目14个，生物过程类的GO条目47个。在细胞组分类别中最显著是“薄膜涂层”，其次是“涂层膜”“黏着斑”条目。在分子功能类别中最显著的是“分子功能”，其次是“连接” “蛋白质结合”条目。在生物学过程类别中最显著的是“生物调节”，其次是“细胞过程的调节”“生物过程的调节”，另外，我们还检测到与脂质合成及代谢相关过程的显著富集，包括“糖异生”“葡萄糖代谢过程”“己糖生物合成过程”“单糖生物合成过程”“糖酵解过程”“丙酮酸代谢过程” “己糖代谢过程”等。 这说明可变剪接对于脂肪沉积调控网络可能有着重要的影响（
图1 GO及KEGG功能富集图
Fig.1 Functional enrichment of GO and KEGG
A：GO富集分析散点图Scatter plot of GO enrichment analysis； B：KEGG富集分析散点图Scatter plot of KEGG enrichment analysis.
对检测到可变剪接事件的基因进行KEGG富集分析，以期研究可变剪接对基因之间互作关系的影响，结果显示有27条通路显著富集，其中直接与脂质合成及代谢相关的有15条通路，例如“脂肪细胞中脂解的调节” “脂肪消化吸收” “脂肪细胞因子信号通路” “PPAR信号通路”“AMPK信号通路” “PI3K-Akt 信号通路”等。这说明可变剪接的参与导致脂肪沉积调控网络的复杂化。
为了获得更直观的注释结果，使用NR、KOG、GO、KEGG和Swissport 5个数据库进行注释，在全部80 756条序列中有69 259条被注释，如
项目 Item | 数据库 Database | 总基因数 Total gene number | 总注释 General notes | |||||
|---|---|---|---|---|---|---|---|---|
| KOG | KEGG | NR | Swiss Port | GO | ||||
| 基因数量 Gene number | 41 004 | 45 984 | 65 902 | 57 880 | 26 653 | 80 756 | 69 256 | |
| 注释比例/% Ratio of notes | 0.78 | 56.94 | 81.61 | 71.67 | 33.00 | - | 85.76 | |
图2 功能注释韦恩图
Fig.2 Venn diagram of functional annotation
图3 NR分类注释统计图
Fig.3 NR classification annotation statistics
KOG数据库注释到41 004条序列在基因区域，占全部序列的50.78%。根据序列信息可以分为26组，其中一般功能预测条目注释到8 406条序列，信号转导机制条目7 221条，注释最少的为未知蛋白组（仅为52条）。从分组可以发现，注释到细胞代谢发育等通用途径的序列占绝大多数。例如：胞内运输、分泌和囊泡运输、脂质运输和代谢、能量生产和转化分别注释到2 822、1 620、920条序列，这些途径证明了在夷陵黄牛脂肪组织中进行着高频率的能量代谢以及物质跨膜运输，这可能与脂肪细胞中脂质积累相关（
图4 KOG 分类注释统计图
Fig.4 KOG classification annotation statistics
GO数据库注释到的基因有26 653个，占全部基因数量的33%，富集到GO条目共1 280个，在分子功能条目上富集到559个条目，30 672条序列，其中以蛋白质结合、核酸结合条目聚集到最多的序列，分别为5 123条和2 706条；生物过程富集到507个条目 15 039条序列，其中转录调控及蛋白磷酸化聚集序列最多，分别为1 495条和1 423条；细胞组分富集到214个条目 7 907条序列，其中膜和膜的组成部分条目聚集的序列最多，分别为1 142条和1 721条。其中本研究检测到与细胞能量储存与代谢等相关条目的显著富集，可能与脂肪组织复杂的生物学过程及活跃的细胞代谢有关（
图5 GO分类注释统计图
Fig. 5 GO classification annotation statistics
KEGG的功能注释结果发现，夷陵黄牛脂肪组织全长转录组数据中共注释到基因45 984个，占全部基因的56.94%，KEGG将注释得到的基因信息分为5个大类：细胞过程、环境信息处理、遗传信息处理、代谢、有机系统，这些类别分别富集到10 745、9 516、5 789、7 005、1 655个基因，如
图6 KEGG 分类注释统计图
Fig.6 KEGG classification annotation statistics
将80 756条序列与Swiss-Prot数据库蛋白库进行比对，共注释得到69 259个基因，占全部基因数量的85.76%。采用 angel 软件，对序列做 CDS 和 Protein 预测，共发现94 458条CDS序列及其对应的pep序列，并且比对到UTR非编码区域的155 214条序列。将数据与基因注释结果对比优化后，得到3' or 5'延长区域共 29 448 处，分别对应21 145条reads，7 176 个基因。
可变剪接已经被证实在95%以上的人类基因中存在，可变剪接的存在使得基因的功能存在不确定性，脂肪调控网络变得更为复
本研究共发现80 756个新基因序列，KOG、KEGG、NR、Swiss Prot、GO数据库分别注释到其中41 004、45 984、65 902、57 880、26 653条序列，在物种同源性比对中发现夷陵黄牛主要与普通牛和瘤牛同源，属于南方黄牛的血统组成，与夏小婷
参考文献References
孙亚红，李俊霞，胥福勋，等.影响肉牛肌内脂肪沉积的遗传、管理和营养因素［J］.中国饲料，2020（2）：102-106.SUN Y H，LI J X，XU F X，et al.Factors related to genetic，management and nutrition affecting fat deposition in beef muscle［J］.China feed，2020（2）：102-106（in Chinese with English abstract）. [百度学术]
万发春，张幸开，张丽萍，等.牛肉品质评定的主要指标［J］.中国畜牧兽医，2004，31（12）：17-19.WAN F C，ZHANG X K，ZHANG L P，et al.Main indexes of beef quality evaluation［J］.Animal science abroad，2004，31（12）：17-19（in Chinese）. [百度学术]
李军涛.日粮VA添加剂量对秦安杂交阉牛牛肉品质及血清生化指标的影响［D］.郑州：河南农业大学，2012.LI J T.Effects of vitamin A on beef quality and serum bioche-mical indicators of Angus × Qinchuan crossbred steers［D］.Zhengzhou：Henan Agricultural University，2012（in Chinese with English abstract）. [百度学术]
UTRERA A R，VAN VLECK L D.Heritability estimates for carcass traits of cattle：a review［J］.Genetics and molecular research：GMR，2004，3（3）：380-394. [百度学术]
MODREK B，LEE C.A genomic view of alternative splicing［J］.Nature genetics，2002，30（1）：13-19. [百度学术]
GALLEGO-PAEZ L M，BORDONE M C，LEOTE A C，et al.Alternative splicing：the pledge，the turn，and the prestige［J］.Human genetics，2017，136（9）：1015-1042. [百度学术]
BLENCOWE B J.The relationship between alternative splicing and proteomic complexity［J］.Trends in biochemical sciences，2017，42（6）：407-408. [百度学术]
YANG X P，COULOMBE-HUNTINGTON J，KANG S L，et al.Widespread expansion of protein interaction capabilities by alternative splicing［J］.Cell，2016，164（4）：805-817. [百度学术]
ZHOU Y，SUN J J，LI C J，et al.Characterization of transcriptional complexity during adipose tissue development in bovines of different ages and sexes［J/OL］.PLoS One，2014，9（7）：e101261［2022-04-18］.https：//doi.org/10.1371/journal.pone.0101261. [百度学术]
VERNIA S，EDWARDS Y J，HAN M S，et al.An alternative splicing program promotes adipose tissue thermogenesis［J/OL］.eLife，2016，5：e17672［2022-04-18］.https：//doi.org/10.7554/eLife.17672. [百度学术]
KAMINSKA D，KUULASMAA T，VENESMAA S，et al.Adipose tissue TCF7L2 splicing is regulated by weight loss and associates with glucose and fatty acid metabolism［J］.Diabetes，2012，61（11）：2807-2813. [百度学术]
CHAO Y H，YUAN J B，LI S F，et al.Analysis of transcripts and splice isoforms in red clover （Trifolium pratense L.） by single-molecule long-read sequencing［J/OL］.BMC plant biology，2018，18（1）：300［2022-04-18］.https：//pubmed.ncbi.nlm.nih.gov/30477428/.DOI： 10.1186/s12870-018-1534-8. [百度学术]
SHIN S C，AHN D H，KIM S J，et al.Advantages of single-molecule real-time sequencing in high-GC content genomes［J/OL］.PLoS One，2013，8（7）：e68824［2022-04-18］.https：//doi.org/10.1371/journal.pone.0068824. [百度学术]
章寒，吴宜静，张一鸣，等．尾叶紫薇全长转录组测序及微卫星标记开发［J］．江西农业大学学报，2021，43（4）：919-930．ZHANG H，WU Y J，ZHANG Y M，et al．Full-length transcriptome sequencing and microsatellite marker development in Lagerstroemia caudata［J］．Acta Agriculturae Universitatis jiangxiensis，2021，43（4）：919-930（in Chinese with English abstract）. [百度学术]
夏丽飞，孙云南，宋维希，等.基于PacBio平台的紫娟茶树全长转录组分析［J］.基因组学与应用生物学，2020，39（6）：2646-2658.XIA L F，SUN Y N，SONG W X，et al.Full-length transcriptome analysis of Zijuan tea（Camellia sinensis var.asssamica（Masters） Kitamura） base on PacBio platform［J］.Genomics and applied biology，2020，39（6）：2646-2658（in Chinese with English abstract）. [百度学术]
张守信，陈英转，郑奋，等.基于PacBio平台的掌叶鱼黄草全长转录组分析［J］.分子植物育种，2022，20（9）：2939-2948.ZHANG S X，CHEN Y Z，ZHENG F，et al.Full-length transcriptome analysis of Merremia vitifolia based on PacBio platform［J］.Molecular plant breeding，2022，20（9）：2939-2948（in Chinese with English abstract）. [百度学术]
SUN M，HUANG D J，ZHANG A L，et al.Transcriptome analysis of heat stress and drought stress in pearl millet based on Pacbio full-length transcriptome sequencing［J/OL］.BMC plant biology，2020，20（1）：323［2022-04-18］.https：//pubmed.ncbi.nlm.nih.gov/32640987/.DOI： 10.1186/s12870-020-02530-0. [百度学术]
ZHANG G Q，SUN M，WANG J F，et al.PacBio full-length cDNA sequencing integrated with RNA-seq reads drastically improves the discovery of splicing transcripts in rice［J］.The plant journal：for cell and molecular biology，2019，97（2）：296-305. [百度学术]
LIU J D，TAN S X，HUANG S Q，et al.ASlive：a database for alternative splicing atlas in livestock animals［J/OL］.BMC genomics，2020，21（1）：97［2022-04-18］https：//doi.org/10.1186/s12864-020-6472-9. [百度学术]
MAO Z，FENG M J，LI Z R，et al.ETV5 regulates hepatic fatty acid metabolism through PPAR signaling pathway［J］.Diabetes，2021，70（1）：214-226. [百度学术]
ZHAO Y，HU X B，LIU Y J，et al.ROS signaling under metabolic stress：cross-talk between AMPK and AKT pathway［J/OL］.Molecular cancer，2017，16（1）：79［2022-04-18］.https：//pubmed.ncbi.nlm.nih.gov/28407774/.DOI： 10.1186/s12943-017-0648-1. [百度学术]
HUANG X J，LIU G H，GUO J，et al.The PI3K/AKT pathway in obesity and type 2 diabetes［J］.International journal of biological sciences，2018，14（11）：1483-1496. [百度学术]
BARTEL D P.microRNAs： genomics，biogenesis，mechanism，and function ［J］.Cell，2004，116（2）：281-297. [百度学术]
MAKDISSY N，HADDAD K，MOUAWAD C，et al.Regulation of SREBPs by sphingomyelin in adipocytes via a caveolin and ras-ERK-MAPK-CREB signaling pathway［J/OL］.PLoS One，2015，10（7）：e0133181［2022-04-18］.https：//doi.org/10.1371/journal.pone.0133181. [百度学术]
YEUNG F，RAMÍREZ C M，MATEOS-GOMEZ P A，et al.Nontelomeric role for Rap1 in regulating metabolism and protecting against obesity［J］.Cell reports，2013，3（6）：1847-1856. [百度学术]
MARIMAN E C M，WANG P.Adipocyte extracellular matrix composition，dynamics and role in obesity［J］.Cellular and molecular life sciences，2010，67（8）：1277-1292. [百度学术]
AI X，HOU X M，GUO T T.C-type natriuretic peptide promotes adipogenic differentiation of goat adipose-derived stem cells via cGMP/PKG/p38 MAPK signal pathway［J］.In vitro cellular & developmental biology-animal，2021，57（9）：865-877. [百度学术]
GUO M，CHEN F，LIN T F，et al.Apelin-13 decreases lipid storage in hypertrophic adipocytes in vitro through the upregulation of AQP7 expression by the PI3K signaling pathway［J］.Medical science monitor：international medical journal of experimental and clinical research，2014，20：1345-1352. [百度学术]
夏小婷，贾媛，赵明，等.夷陵黄牛mtDNAD-loop区遗传多样性研究［J］.中国牛业科学，2017，43（1）：4-7.XIA X T，JIA Y，ZHAO M，et al.Genetic diversity of mitochondrial DNA D-oop region in Yiling cattle［J］.China cattle science，2017，43（1）：4-7（in Chinese with English abstract）. [百度学术]
YANG J，WARISS H M，TAO L D，et al.De novo genome assembly of the endangered Acer yangbiense，a plant species with extremely small populations endemic to Yunnan Province，China［J/OL］.GigaScience，2019，8（7）：giz085［2022-04-18］.https：//doi.org/10.1093/gigascience/giz085. [百度学术]
BYRNE A，COLE C，VOLDEN R，et al.Realizing the potential of full-length transcriptome sequencing［J/OL］.Philosophical transactions of the Royal Society of London.Series B，Biological sciences，2019，374（1786）：20190097［2022-04-18］.https：//doi.org/10.1098/rstb.2019.0097. [百度学术]