多浪羊与小尾寒羊皮下脂肪组织转录组分析

刘天义，冯 卉，Salsabeel Yousuf，解领丽，苗向阳

(中国农业科学院北京畜牧兽医研究所，北京 100193)

绵羊是全世界主要的肉、奶、皮毛的牲畜资源。羊肉是我国日常膳食结构的重要组成部分。随着人们生活水平的提高，人们对于羊肉品质的要求也随之提高，而脂肪含量对肉品质具有重要影响[1-2]。研究脂肪组织的脂肪沉积机制不仅可以为人们提供高品质肉产品，还可以提高养殖业饲料转化率，避免浪费。同时通过对脂肪组织的研究还能为由肥胖引起的相关疾病的预防和治疗提供一定依据[3-5]。

脂肪组织遍布哺乳动物全身，主要分为白色脂肪组织和棕色脂肪组织两类[6-7]。研究表明，脂肪合成与分解代谢过程受到C/EBPs、PPARγ、BMP4等转录因子[8-10]，以及PPAR、MAPK、TGFβ、Insulin等信号通路的调控[11-12]。Xu等[13]研究发现，SCD、FAS、FAP4等基因在脂肪酸代谢过程中发挥作用，其中，SCD1是SCD的一种亚型，可作为脂肪沉积的潜在标志物，SCD1是SREBP1调控的下游基因，当SCD1的启动子上的固醇调节原件与SREBP1c转录因子结合后上调转录，脂质合成增加[14]。RNA是遗传信息的重要载体，mRNA可以编码蛋白实现遗传信息的传递，lncRNA和miRNA均可以靶向mRNA来调控基因的表达[15-17]。转录组测序技术由于具有高灵敏度，已经被广泛应用。Miao等[18]通过转录组测序手段，对多塞特羊与小尾寒羊脂肪组织进行转录组测序，以更好地了解两品种绵羊脂肪代谢的生物学机制。Huang等[19]通过对和牛与荷斯坦牛皮下脂肪组织进行转录组测序，筛选得到差异表达基因，进一步对两品种牛脂肪沉积的分子机制进行解析。为更深入了解mRNA在绵羊脂肪组织中的作用，本研究选用脂肪沉积存在差异的多浪羊与小尾寒羊皮下脂肪组织为试验材料，其中多浪羊脂臀较大，属于肉脂兼用型绵羊，小尾寒羊属于短瘦尾型绵羊，肉用性能较差，皮下脂肪较少，而脂肪主要分布在内脏周围，以此来适应恶劣的环境条件，二者为研究脂肪的沉积机制提供了良好的试验材料。目前转录组测序技术较为成熟，以往虽有通过转录组技术来对脂肪组织进行的研究报道，但对于多浪羊与小尾寒羊皮下脂肪组织沉积的生物学差异还有待于深入研究。

本研究通过对多浪羊与小尾寒羊皮下脂肪组织进行转录组测序分析，得到两品种羊皮下脂肪组织中差异表达的mRNA，并从中筛选出与脂肪沉积相关的候选基因。

1 材料与方法

1.1 试验动物和样品制备

本试验以脂肪沉积存在差异的多浪羊和小尾寒羊为研究对象，选取健康无病、体况良好，种内个体体重相近(约50 kg)的雌性成年多浪羊和小尾寒羊各3只。采集背最长肌皮下脂肪组织为试验材料，迅速将皮下脂肪组织放入液氮冻存。该过程使用已消毒器械，减少RNA的降解。试验设置多浪羊皮下脂肪组织和小尾寒羊皮下脂肪组织两组样本，每组样本3个重复，分别为D-PF-1、D-PF-2、D-PF-3、X-PF-1、X-PF-2、X-PF-3。

1.2 方法

1.2.1 脂肪组织RNA的提取及质控 取适量的脂肪组织，采用TRIzol(Invitrogen)法提取脂肪组织中的总RNA。并使用DNase I (RNase-free，TaKaRa)去除DNA。利用1%琼脂糖凝胶电泳与ND-2000(NanoDrop Technologies)对所提RNA的浓度、纯度以及完整性进行初步检测，再用2100 Bioanalyzer(Agilent)进行检测，记录OD260 nm/OD280 nm以及RIN值。

1.2.2 脂肪组织cDNA文库构建和RNA测序 将质检合格的总RNA进行文库构建。采用TruSeqTMstranded total RNA Kit(Illumina，San Diego，CA)试剂盒。首先从2 μg总RNA中除去rRNA；加入Fragmentation Buffer将RNA片段化；采用SuperScript double-stranded cDNA synthesis kit (Invitrogen,CA)试剂盒，加入六碱基随机引物(Illumina)，利用随机引物，以mRNA为模板反转录合成第一链cDNA，随后合成第二链cDNA，形成稳定的双链结构；加入End Repair Mix将粘性末端补成平末端，随后在3′末端加上一个A碱基，用于Y字形的接头，扩增前用UNG酶消化cDNA第二链；cDNA 经过 PCR 富集(sample preparation Kit(Illumina，San Diego,CA)后，DNA clean beads筛选 200～300 bp 的条 带。经 TBS380(Picogreen)定量后，文库使用 Illumina NovaSeq 6000 测序平台进行高通量测序。

1.2.3 脂肪组织mRNA测序分析 为保证后续分析的准确性，需要使用fastp[20]软件(https://github.com/OpenGene/fastp)将原始数据(raw data)中包含的测序接头序列、低质量读段、不确定碱基信息率较高序列及长度过短序列进行过滤，从而得到高质量的测序数据(clean data)。随后使用Hisat2[21](https://daehwankimlab.github.io/hisat2/)将筛选后的数据定位到参考基因组，获取用于后续分析的有效数据，同时将对比的结果用分析软件RSeQC-2.3.6[22-23]进行质量评估。在已有的参考基因组基础上使用软件stringTie[24-25](https://ccb.jhu.edu/software/stringtie/)对每个样本进行组装拼接，与已有的转录本进行比较，获得没有注释的新转录本。用软件RSEM[26-27](http://deweylab.github.io/RSEM/)进行mRNA的表达量计算，再用软件DESeq2筛选样本之间差异表达的mRNA，将条件设定为|Fold change|≥2，Padjust≤0.05进行筛选。

1.2.4 脂肪组织差异表达mRNA功能富集分析 应用Cytoscape中的Clue-GO[28]软件将差异表达的mRNA进行GO功能注释和KEGG通路富集分析，GO是基因本体论联合会建立的数据库(geneontology，http://geneontology.org/)，分为细胞组分(cellular component，CC)，分子功能(molecular function，MF)和生物学过程(biological process，BP)。KEGG[29]京都基因和基因组百科全书(https://www.genome.jp/kegg/)是系统分析基因功能、联系基因组信息和功能信息的数据库，可将基因按照参与的pathway通路或行使的功能分类。将集中的基因显示在KEGG通路图上，展示其参与的KEGG注释通路图。利用多重检验校正方法BH(Benjamini and Hochberg)对P值进行校正，当P值(Padjust)小于等于0.05时则显著富集。

1.2.5 脂肪组织mRNA蛋白网络互作分析 将GO和KEGG注释得到的信息，结合文献筛选出与脂肪发育相关的条目和信息通路，选择其中涉及到的基因，利用Cytoscape[30]软件制作网络图，找出网络中处于节点位置的比较重要的基因，在文献中找出这些基因具体的功能以及在脂肪组织发育中的作用。

1.2.6 qRT-PCR验证 采用实时荧光定量PCR验证基因的表达水平。利用PerfectStartTMGreen qPCR SuperMix试剂盒在LightCycler480 Ⅱ Real-time PCR Instrument(Roche，Swiss)上配置PCR反应体系：2×PerfectStartTMGreen qPCR SuperMix 5 μL，10 μmol·L-1上、下游引物各 0.2 μL，cDNA 1 μL，Nuclease-free H2O 3.6 μL。PCR程序：94 ℃ 30 s；94 ℃ 5 s，60 ℃ 30 s，45个循环。循环结束后利用熔解曲线检测产物特异性：从60 ℃缓慢升温至97 ℃，每1 ℃采集5次荧光信号。利用2-ΔΔCt法进行样本间表达量的计算[31]，用t-检验统计分析相对表达量，数据结果表示为“平均数±标准差(Mean±SD)”，P<0.05为差异显著。对应的引物序列见表1。

表1 基因以及对应的引物Table 1 Gene and the corresponding primer sequences

2 结 果

2.1 多浪羊与小尾寒羊皮下脂肪组织测序数据分析

将原始测序数据总条目数(raw reads)中质检不合格的条目去除后得到质控后测序数据的总条目数(clean reads)，多浪羊与小尾寒羊皮下脂肪组织clean reads分别为135 393 598、140 587 412、142 895 188、136 746 686、128 433 640和130 999 204。将这些reads与参考基因组对比得到有效数据(mapped reads),以及在参考序列上有多个比对位置的clean reads数目(multiple mapped)和唯一比对位置的clean reads数目(unique mapped)见表2。Mapping率高于94%，说明参考基因组注释完整，且相关试验不存在污染情况，为后续数据分析做了很好的铺垫。

表2 6个样本比对参考基因组结果Table 2 The result of 6 samples after mapping to the reference genome

2.2 已知与新mRNA分析及表达量计算

将数据进行已知和新mRNA筛选，结果显示，在多浪羊和小尾寒羊皮下脂肪组织样本中已知和新mRNA的数量分别为38 672和1 606个。每组设置3个生物学重复计算mRNA的表达量，选取|Fold change|≥2，Padjust≤0.05作为筛选差异表达mRNA的阈值。发现两组差异表达的mRNA共839个(已知mRNA 828个，新mRNA 11个)，如图1所示，其中在多浪羊组中上调和下调的差异表达基因数目分别为320和519个。

图1 差异表达基因火山图Fig.1 Volcano map of differentially expressed genes

2.3 脂肪组织差异表达mRNA的功能富集分析

对差异表达的mRNA进行GO和KEGG分析。GO注释结果显示，在生物学过程中，差异表达的mRNA主要富集于脂质分解代谢过程、脂质生物合成过程、脂质分解代谢负调控过程、MAPK级联反应调控、低密度脂蛋白颗粒介导的信号传导、对甘油三酯的反应、脂质代谢过程、MAPK级联反应负调控等条目；在分子功能方面，主要富集于细胞外基质结构成分、心磷脂结合、甘油酸激酶活性、磷脂酰甘油结合等条目；在细胞组分，主要富集于细胞外基质、脂蛋白颗粒、血浆脂蛋白颗粒等条目(图2)。在KEGG通路富集中也富集到与脂代谢相关的通路，如PI3K-Akt、MAPK、胰岛素以及PPAR等信号通路(图3)。

图2 差异表达基因GO富集分析Fig.2 GO enrichment analysis of the differentially expressed genes

图3 差异表达基因KEGG通路富集Fig.3 KEGG pathway enrichment analysis of the differentially expressed genes

2.4 脂肪组织差异表达基因蛋白互作网络图的构建

构建与脂质代谢相关差异表达基因蛋白互作网络，能够更直观的展示多浪羊和小尾寒羊与脂肪沉积相关的基因间的调控关系(图4)。本研究以网络中各个节点的联通性(Degree)为依据，筛选出关键节点[32-33]。Degree值越大，节点越大，表明该节点越关键。在蛋白互作网络中发现，COL1A1、AKT3、AKT2、COL1A2、COL3A1等基因编码的蛋白处于关键节点位置，可能在皮下脂肪代谢和成脂分化中具有重要调控作用。

2.5 qRT-PCR验证

为了验证测序结果，随机选择6个差异表达的mRNAs进行qRT-PCR验证，结果如图5所示，MGST3、COL1A1、AKT3在多浪羊皮下脂肪组织中下调，PCK1、PPP2R5A、LOC101113583在多浪羊皮下脂肪组织中显著上调。以上结果与测序结果一致，表明测序结果可靠。

图5 差异表达基因qRT-PCR验证Fig.5 qRT-PCR validation for differentially expressed genes

3 讨 论

mRNA是指导蛋白质合成的模板，在生物学过程中发挥着重要作用[34]，可为培育优质畜禽品种以及研究和预防与脂代谢相关的疾病提供依据。本研究对多浪羊与小尾寒羊皮下脂肪组织中的mRNA进行测序和分析，共鉴定到38 672个mRNAs，其中差异表达的mRNAs共839个，在多浪羊中有320个上调表达，519个下调表达。通过qRT-PCR对MGST3、COL1A1、AKT3、PCK1、PPP2R5A、LOC101113583进行了验证，与测序结果一致，表明测序结果可靠。

有研究表明，细胞外基质可以影响脂肪的形成，而胶原是脂肪细胞外基质的主要成分，对脂肪的形成具有一定的调节作用[35]，推测研究编码胶原的基因COL1A1对于脂肪的沉积具有重要意义。在秦川牛中的研究表明，COL1A1的表达量与肉质存在一定的相关性[36]。COL1A1在高脂饮食的肥胖人皮下脂肪组织中和高脂饲喂小鼠肌腱内的表达量均降低[37-38]，与本研究趋势一致，表明COL1A1可能参与了脂肪代谢的调控。本研究中，AKT2显著富集于PI3K/Akt信号通路中。PI3K/Akt是胰岛素信号传导的主要通路，通过调控胰岛素水平参与糖脂代谢，可以抑制脂类分解。AKT2是胰岛素信号通路的关键激酶，仅在胰岛素靶组织中高表达，包括脂肪、肝、肌肉等组织，而胰岛素可以调节脂肪细胞的分化，即推测AKT2基因可能调控了脂肪细胞的分化。敲除AKT2基因的小鼠表现出脂肪组织减少、胰岛素抵抗和轻度糖尿病等现象[39]，本研究中，AKT2在多浪羊皮下脂肪组织中高表达，表明AKT2基因对多浪羊脂类细胞分解具有促进作用。综上所述，COL1A1与AKT2不仅在脂肪沉积过程中起到了调控作用，而且还可能成为治疗糖尿病等相关代谢疾病的潜在靶基因。

节点的大小代表连通性大小，白色表示下调，灰色表示上调The size of the node represents the level of Degree,white means down-regulation and gray means up-regulation图4 差异表达基因蛋白互作网络图Fig.4 Network diagram of differentially expressed gene proteins interaction

通过差异表达基因GO注释和KEGG富集分析表明，差异表达基因主要参与了PPAR信号通路、MAPK信号通路、脂肪细胞因子信号通路、脂质代谢过程等与脂类代谢相关的生物学过程。本研究鉴定到SCD、LPL、PCK1、PPP2R5A等基因与脂类代谢、脂类生物合成过程以及胰岛素信号通路、PPAR信号通路过程相关。PPAR信号通路是与脂肪酸的生物合成和代谢以及成脂分化相关的通路[40-41]，本研究中，PPAR信号通路显著富集了10个差异表达基因，其中SCD、LPL和PCK1在多浪羊皮下脂肪组织中上调表达，SCD是一种内质网结合酶，是参与单不饱和脂肪酸生物合成的关键酶，并且具有多种亚型。由于PPARγ2对于甘油三酯合成以及脂肪酸吸收的相关基因具有调控作用，所以在脂质代谢过程中具有重要作用，SCD1可以作为其靶基因在脂代谢过程中发挥作用。SCD1的表达水平受类固醇激素、甲状腺激素、共轭亚油酸等影响。LPL是脂肪分解的关键限速酶，有利于脂肪的沉积[42-45]，并且主要在转录后蛋白水平起作用。LPL作为脂肪细胞的特征蛋白质之一，在脂肪细胞分化过程中具有重要作用，可以充当前体脂肪细胞向成熟脂肪细胞分化过程的标志。本研究中，LPL在多浪羊皮下脂肪组织中上调表达，与前人研究结果一致，推测LPL同样对绵羊脂肪的生成具有促进作用。PCK1是调节糖异生和甘油生成的主要酶之一[46]，与肥胖相关，能通过GTP催化草酰乙酸(oxaloacetic acid，OAA)生成磷酸烯醇式丙酮酸(phosphoenolpyruvate，PEP)、GDP以及二氧化碳，在甘油生成以及糖异生过程中均可见到。在白色脂肪组织中，PCK1可以调节游离脂肪酸酯化生成甘油三酯，PPARγ是其启动子中的结合位点，与野生型小鼠相比，PPAR缺失小鼠患有严重胰岛素抵抗[47]，还有研究表明，敲除小鼠PCK1基因会导致糖尿病以及肥胖等[48]。因此PCK1对于维持脂质代谢至关重要。

PPP2R5A是蛋白磷酸酶2A的一个亚单位，与细胞增殖分化相关,可以通过调节CDK1/2和CHK2的活性促进有丝分裂，并且可以调节脂质激酶的活性[49-50]。本研究中，PPP2R5A在多浪羊皮下脂肪组织中上调，GO注释发现其主要参与了细胞过程和脂质代谢过程的调节。PPP2R5A显著富集于PI3K-Akt信号通路，这是细胞内重要的通路，可以调控细胞生长增殖、分化等过程，因此推测PPP2R5A可能在脂肪细胞增殖分化等过程中发挥作用。MGST3与LOC101113583在脂肪沉积的研究中较为少见，但在本研究中显著富集到脂质生物合成与代谢生物学过程中，对其功能的研究还需进一步验证。

本研究对多浪羊与小尾寒羊皮下脂肪组织中差异表达的mRNA进行了分析，结果表明，这些差异表达基因可能参与脂代谢过程，并且从中筛选到影响脂肪沉积的候选基因以及相关通路。但对于筛选到的候选mRNA在生物学功能以及调控的分子网络还需要进一步探索。

4 结 论

本研究鉴定分析了多浪羊和小尾寒羊皮下脂肪组织中mRNA的表达情况，发现差异表达的mRNAs共839个，在多浪羊皮下脂肪组织中320个上调，519个下调；通过对差异表达基因进行GO和KEGG富集分析，筛选出COL1A1、AKT2、SCD、LPL、PCK1与PPP2R5A等与绵羊脂质代谢相关的基因，这些结果可为进一步理解多浪羊与小尾寒羊皮下脂肪沉积与脂代谢差异机制提供理论基础。