动物脂肪细胞分化过程中全转录组3´-UTR动态变化研究

2022-03-30杜志强赵素娟陈耀峰殷宗俊杨彩侠

东北农业大学学报 2022年3期

杜志强，赵素娟，陈耀峰，丁然，殷宗俊，杨彩侠

（1.东北农业大学动物科学技术学院，哈尔滨 150030；2.长江大学动物科学学院，湖北荆州 434025;3.安徽农业大学动物科学技术学院，合肥 230036）

动物脂肪组织主要由脂肪细胞构成，通过储存能量和执行内分泌功能影响新陈代谢，参与调控动物健康和肉品品质[1]。脂肪生成关键基因和转录因子构成的基因网络调控脂肪细胞分化过程，其中过氧化物酶体增殖物激活受体γ（PPARγ）和CCAAT增强子结合蛋白α（C/EBPα）是最主要的两类调控因子[2]。不同部位脂肪组织生理生化特性不同，糖利用、脂代谢机制和激素敏感性也存在差异，如肌内脂肪含量是影响肉质重要因素，腹部脂肪含量影响屠宰效率、皮下白色脂肪组织，发挥储存脂肪和抗炎作用等[3]。

比较分析发现，哺乳动物（人和小鼠）脂肪细胞系和家禽前脂肪细胞分化存在异同。与哺乳动物相比，禽类脂肪细胞分化同样受到PPARγ与C/EBPα等因子影响，但其表达模式及调控作用均有差异。C/EBPα和PPARγ发挥启动禽类前脂肪细胞分化的关键作用，且C/EBPα调控PPARγ转录活性，而PPARγ对C/EBPα表达无直接调控作用[4]。脂肪型脂肪酸结合蛋白（A-FABP）一般在哺乳动物脂肪细胞分化中晚期高表达，在鸡脂肪细胞分化早期表达量明显上升，比C/EBPα等基因更早表达[5]。近年来，在深入研究和比较不同动物脂肪生成分子机制上，转录组测序分析发现众多关键转录因子和信号通路[6-9]。

真核生物有多种mRNA转录后调控机制，其中mRNA 3´-非翻译区（3´-UTR）参与发挥重要作用。前体mRNA成熟过程中，3´-UTR发生剪切和多聚腺苷酸化反应，而不同多聚腺苷酸化位点选择，导致转录异构体具有不同3´-UTR长度和序列组成[10]。3´-UTR长度变化导致丢失或获得顺式作用元件，如富含AU序列元件（ARE），富含GU序列元件（GRE）等，可被microRNA、RNA结合蛋白（RBP）等反式作用因子结合，进而调控mRNA亚细胞定位、稳定性及翻译效率[11]。因此，3´-UTR具有以下主要功能：通过顺式作用元件调节mRNA稳定性；调控mRNA定位表达和翻译；结合多个RBP行使调控功能；调控蛋白-蛋白相互作用[11-12]。

随3´-UTR功能研究深入，发现其与动物生产性能关系密切。首先，运用候选基因法，即挑选功能已知相关基因，通过生物信息学分析、关联性分析和构建双荧光素酶报告基因表达系统等方法，检测3´-UTR同性状关联。Bartz等发现ME1基因3´-UTR中有7个新变异，通过关联分析发现其多态性与波兰长白猪和大白猪背膘厚度和平均日增重显著相关[13]。其次，通过分析microRNA与其靶mRNA 3´-UTR序列互补配对，可探究3´-UTR如何影响基因表达。Li等发现，miRNA-122可靶向于Vanin 1基因3´-UTR，调控Vanin 1表达，影响鸡肝脏脂质代谢过程[14]。但前脂肪细胞分化过程中是否存在3´-UTR动态变化，目前还未见相关报道。

本研究通过使用Dapars和APAtrap两个选择性多聚腺苷酸化（APA）预测软件，分析四种动物（鸡、鸭、猪和小鼠）不同分化时期前脂肪细胞RNA-seq数据，探究脂肪生成过程中转录组水平上3´-UTR长度存在动态变化规律。研究将有助于理解动物前脂肪细胞分化过程中转录调控机制，为动物脂肪生长发育生物学机制研究提供新视角。

1 材料与方法

1.1 转录组测序数据

本研究使用4个不同动物前脂肪细胞不同分化时期4套转录组测序（RNA-seq）数据集，包括京海黄鸡肌肉内脂肪、北京鸭皮下脂肪和长白猪皮下脂肪来源的前脂肪细胞，以及小鼠3T3-L1细胞[7-9]。数据均从国家生物技术信息中心（NCBI）短序列档案（SRA）数据库中下载，SRA号分别为SRX2037700（鸡）、SRX4646736（鸭）、SRX5121442（猪）和SRX345947（鼠）。京海黄鸡有12个样本，分别在前脂肪细胞分化0、2、4、6 d进行测序（记为I0、I2、I4和I6），每个时间点3个重复；北京鸭36个样本分别在前脂肪细胞分化-48、0、12、24、48、72 h进行测序(记为T-48、T0、T12、T24、T48和T72），每个时间点6个重复；长白猪共12个样本，分别在前脂肪细胞分化0、2、4、8天进行测序（记为I0、I2、I4和I8），每个时间点3个重复；3T3-L1细胞测序8个样本，时间点分别为3T3-L1细胞指数扩散期、48 h融合期、诱导分化10 h和第6天（记为EP、C、T10 h和T6 d），每个时间点2个重复。

1.2 生物信息学分析

4套RNA-seq数据通过fastqc网址控制质量。将clean reads通过HISAT2[7]比对到参考基因组（Gallus_gallus.GRCg6a（鸡）；Anas_platyrhynchos_platyrhynchos.CAU_duck1.0（鸭）；Sus_scrofa.Sscrofa 11.1（猪）；GCF_000001635.27_GRCm39_genomic.fna（小鼠）。最后，使用samtools将sam文件处理成bam文件，使用featureCounts计算基因的表达量。

1.3 全转录组水平3´-UTR长度差异分析

Dapars（Dynamic Analyses of Alternative PolyAdenylation from RNA-Seq，version 0.9.1)通过比较RNA-seq数据推断APA位点[15]。Dapars将注释的3´-UTR末端作为远端poly（A）位点，据此预测其上游近端poly（A）位点，推断APA差异使用。因此，通过比较RNA-seq数据，Dapars可推断出APA起始位置，长或短3´-UTR转录本表达水平，从而确定APA动态用法。转录本应用APA的情况，可通过比较两组数据之间PDUI变化预测：PDUI_Group_diff值为正，应用近端poly（A）位点（对应短3´-UTR）；PDUI_Group_diff值为负，应用远端poly（A）位点（对应长3´-UTR）。

APAtrap也可从RNA-seq数据鉴定和量化poly（A）位点，鉴定不同条件下具有显著APA动态使用的基因[16]。APAtrap识别新的或延长3´-UTR，并精确预测所有潜在poly（A）位点。APAtrap可通过整合一个基因中所有poly（A）位点位置和表达水平，探索不同使用模式APA事件。通过比较两组数据之间r变化，预测转录本应用APA的情况：r值为正，则应用近端poly（A）位点（对应短3´-UTR）；r值为负，则应用远端poly（A）位点（对应长3´-UTR）。

使用Dapars分析时，PDUI参数设为0.25，其他参数为默认参数，按照adjusted.P_val＜0.05和Pass_Filter=Y进行筛选；使用APAtrap分析时，所有参数为默认参数，按照adjusted P-value＜0.05和PD≥0.2筛选。

1.4 基因表达水平和3´-UTR长度差异分析

1.5 功能富集分析

对表达水平和3´-UTR长度有显著差异基因（DE-exp基因和DE-3´-UTR基因）进行Gene Ontology（GO）功能注释和Kyoto Encyclopedia of Genes and Genomes（KEGG）通路富集分析。GO富集结果包括3部分：分子功能（Molecular function，MF）、生物过程（Biological process,BP）和细胞组成（Cellular component，CC）。每部分富集结果（P＜0.05）仅在图中展示前5个（不足5个的全部展示），并列出富集条目较多的生物过程和与脂肪生成相关的生物过程。此外，P＜0.05的生物功能和信号通路被认为显著富集，并绘制成KEGG气泡图。

1.6 microRNA调控分析

3´-UTR长短不同，导致其序列成分存在差异，影响其功能调控，如microRNAs同3´-UTR结合与否，通常影响翻译抑制、蛋白表达降低和/或转录本降解[17]。通过TargetScan（http://www.targetscan.org/vert_72/）和miRbase数据库（http://www.mirbase.org/search.shtml）确定和分析靶向DE-3´-UTR基因microRNA。

2 结果与分析

2.1 RNA-seq数据初步分析

4种动物RNA-seq数据比对和质量控制结果如下：京海黄鸡12个测序文库，共产生1 382 347 320的原始短序列（rawreads），质控后获得1 328 693 812的clean reads（每个文库比例：95.81%~96.52%）；clean reads比对率，78.62%~92.64%。北京鸭36个测序文库，1 908 767 014 raw reads，质控后获得1 855 198 888的clean reads（每个文库比例：95.96%~97.78%）；clean reads比对率，69.23%~79.72%。长白猪12个测序文库，382 731 382 raw reads，质控后获得375 439 840的clean reads（每个文库比例：97.52%~98.37%）；clean reads比对率，96.40%~96.97%；3T3-L1细胞8个测序文库，共产生1 579 675 811 raw reads，质控后获1 434 061 091的clean reads（每个文库比例：85.91%~95.94%）；clean reads比对率，91.91%-97.42%。

2.2 前脂肪细胞分化过程中3´-UTR的长度动态变化

使用Dapars和APAtrap两个软件分析RNA-seq数据，对前脂肪细胞不同分化时期进行两两比较，寻找有3´-UTR差异基因，并分析其长度动态变化。

2.2.1 京海黄鸡

分析京海黄鸡肌内前脂肪细胞测序数据，Dapars分别得到352、440和111个DE-3´-UTR基因（I0vsI2、I2vsI4、I4vsI6），APAtrap分别得到1 378、1 411和357个DE-3´-UTR基因（I0vsI2、I2vsI4、I4vsI6）。两个软件均预测到共同DE-3´-UTR基因数分别为168、230和35个。

与Dapars相比，APAtrap在识别功能重要的APA事件和不同条件下动态APA使用方面，功能更为强大，准确性更高[16]。进一步分析APAtrap的3 146个DE-3´-UTR基因长短转录本使用情况，发现前脂肪细胞分化后第2天，3´-UTR显著变长；第4天与第2天相比，依然存在3´-UTR变长现象；然而，第6天与第4天相比，3´-UTR长度变短。此外，3´-UTR变短与3´-UTR变长的基因数比值随前脂肪细胞分化不断升高（见图1A）。由此可见，京海黄鸡前脂肪细胞在分化前期（I0-I4）倾向于使用3´-UTR更长的转录本，分化后期倾向于使用3´-UTR更短的转录本，分化过程倾向于使用3´-UTR更短的转录本。比较I2vsI4、I4vsI6的DE-3´-UTR基因表达水平，发现使用短3´-UTR基因表达水平更高（见图1B），可能是短3´-UTR逃避microRNA结合，导致基因表达水平更高。

2.2.2 北京鸭

北京鸭皮下前脂肪细胞RNA-seq数据分析，得到DE-3´-UTR基因数为：Dapars，37、161、16、40和42个基因（T-48vsT0、T0vsT12、T12vsT24、T24vsT48、T48vsT72）；APAtrap，10、426、0、0和0个基因（T-48vsT0、T0vsT12、T12vsT24、T24vsT48、T48vsT72）；共同基因数分别为0和17个。

由于APAtrap在T12vsT24、T24vsT48、T48vsT72各组比较中均未预测到DE-3´-UTR基因，因此增加T0vsT24、T0vsT48、T0vsT72组比较，分别得到564、402和478个DE-3´-UTR基因。分析APAtrap得到的1 880个DE-3´-UTR基因长短转录本使用情况，发现在增殖和分化过程中均倾向于使用更短的转录本，且在前脂肪细胞分化24 h后，3´-UTR变短的基因数最多（见图1C）。由此可见，北京鸭前脂肪细胞分化也倾向于使用更短的转录本。比较DE-3´-UTR基因表达水平，发现在T-48vsT0、T0vsT12和T0vsT72差异基因中，使用短3´-UTR的基因表达水平更高（见图1D）。

2.2.3 长白猪

长白猪皮下前脂肪细胞RNA-seq测序数据分析，Dapars和APAtrap两个软件均未预测到有3´-UTR差异基因。

2.2.4 小鼠3T3-L1

在3T3-L1细胞测序数据中，由Dapars分析得到的DE-3´-UTR基因数为7、25和32个基因（EPvsC、CvsT10 h、T10 hvsT6 d）；由APAtrap分析得到的DE-3´-UTR基因数为124、46和125个基因（EPvsC、CvsT10 h、T10 hvsT6 d）；在各组比较中，两个软件均预测到DE-3´-UTR基因数分别为3、9和18个。

分析APAtrap得到的295个DE-3´-UTR基因长短转录本使用情况，发现3T3-L1细胞在增殖阶段，基因3´-UTR显著变长；但是诱导分化10 h后，基因3´-UTR变短；分化的第6天与10 h相比，基因3´-UTR长度进一步显著变短（见图1E）。由此得出结论，3T3-L1细胞分化过程倾向于使用更短的转录本。比较DE-3´-UTR基因表达水平，发现Cvs10h组，使用长3´-UTR的基因表达水平更高，而10hvs6d组，使用短3´-UTR的基因表达水平更高（见图1F）。

图1 京海黄鸡肌内前脂肪细胞分化过程中3´-UTR长度动态变化Fig.1 Dynamic changesof 3´-UTR length during intramuscular preadipocytedifferentiation in Jinghai Yellow chicken

2.3 差异基因功能富集分析

2.3.1 京海黄鸡

对3 146个DE-3´-UTR基因进行GO和KEGG功能富集分析（个别富集到的条目很少甚至没有，因此未绘制直观图，见图2）。其中，GO生物过程分析发现：I0vsI2显著富集于蛋白质定位调节、脂质代谢与合成和细胞自噬等；I2vsI4主要集中在核苷酸代谢和生物合成、蛋白磷酸化脱磷酸化和转化生长因子受体信号通路；I4vsI6主要集中在基因表达、核苷磷酸生物合成与代谢。此外，KEGG通路分析发现：I2vsI4主要富集在磷脂酰肌醇信号系统、磷酸肌醇代谢和FoxO信号通路。

图2 京海黄鸡肌内前脂肪细胞不同分化时期DE-3´-UTR基因功能富集和通路分析Fig.2 Functional enrichment and pathway analysis of DE-3´-UTR genes at different differentiation stages of intramuscular preadipocytes in Jinghai Yellow chicken

GO功能富集和通路分析DE-exp基因发现，与DE-3´-UTR基因显著富集到的生物过程和通路存在显著不同。GO生物过程：I0vsI2组主要集中在脂质代谢、细胞黏附、细胞周期和细胞骨架组织等；I2vsI4组，除细胞黏附、细胞周期和细胞骨架组织外，还有细胞凋亡、蛋白质改性过程和细胞因子产生等；I4vsI6组，除脂质代谢、细胞周期、细胞骨架组织和细胞凋亡信号外，还有MAPK级联调控、蛋白激酶B信号和脂肪细胞分化调控等。KEGG通路分析：显著富集的主要通路是细胞周期、DNA复制和细胞老化。

根据GO和KEGG分析结果和已有文献报道，筛选得到部分DE-3´-UTR基因作为靶基因用于进一步研究（见表1）。

2.3.2 北京鸭

对北京鸭1 880个DE-3´-UTR基因进行GO和KEGG功能富集分析（见图3）。GO富集（P＜0.05）得到BP和MF条目较少，展示排名最前的5个；KEGG通路富集分析（P＜0.05）得到的个别条目很少甚至没有，不绘制直观图。GO分析：T-48vsT0，主要是核苷酸结合等分子功能；T0vsT12和T0vsT72，主要集中在系统发育和大分子生物合成过程；KEGG通路分析：T0vsT12，主要集中在TGF-beta信号通路、焦点粘连、甘油磷脂新陈代谢和MAPK信号通路；T0vsT72，主要集中在TGF-beta信号通路和磷酸肌醇代谢。

图3 北京鸭皮下前脂肪细胞不同分化时期DE-3´-UTR基因功能富集和通路分析Fig.3 Functional enrichment and pathway analysisonDE-3´-UTR genes at different differentiation stages of subcutaneouspreadipocytedifferentiation in Pekin duck

DE-exp基因GO功能富集和KEGG通路分析发现，T-48vsT0组，主要集中在信号传导和生物调节生物过程；T0vsT12组，主要有核糖体、Wnt信号通路、TGF-beta信号通路和细胞周期等信号通路；T0vsT72组，有焦点粘连、细胞周期和脂肪酸代谢等信号通路。与DE-3´-UTR基因显著富集到的通路相比，仅存在TGF-beta信号通路和焦点粘连两条相同通路。由GO和KEGG分析结果，筛选得到部分DE-3´-UTR基因作为靶基因（见表1）。

表1 靶基因及其功能Table1 Target genesand their functions

2.3.3 3T3-L1细胞

小鼠3T3-L1细胞分化相关的295个DE-3´-UTR基因（见图4），GO功能富集分析后发现显著富集的生物过程包括：EPvsC，MAPK活性负调控、葡萄糖运输等；Cvs10h，主要集中在葡萄糖代谢过程正调控和细胞脂质代谢等生物过程；10hvs6d，主要集中在蛋白质定位、线粒体自噬和细胞周期等。KEGG通路富集分析发现：EPvsC，主要集中在mTOR、EGFR酪氨酸激酶抑制剂和Ras信号通路。

图4 3T3-L1细胞不同分化时期DE-3´-UTR基因功能富集和通路分析Fig.4 Functional enrichment and pathway analysisonDE-3´-UTR genesat different differentiation stagesof 3T3-L1 cells

DE-exp基因GO功能富集和KEGG通路分析发现：EPvsC组，主要集中在甘油脂代谢、细胞黏附和细胞周期等生物过程，以及蛋白质消化吸收通路。Cvs10h组，有脂质定位、蛋白定位调节、神经元投射发育和mRNA加工等生物过程；10hvs6d组，除脂质代谢、细胞周期、神经元投射发育和mRNA加工等生物过程，及细胞自噬和Wnt信号通路等。由此可见，DE-exp基因与DE-3´-UTR基因富集到显著不同的生物过程和信号通路。根据GO和KEGG分析结果，筛选得到部分DE-3´-UTR基因作为靶基因用于进一步研究（见表1）。

2.4 microRNA调控

通过TargetScan和miRbase数据库搜索数个DE-3´-UTR基因microRNA结合位点（具有较多microRNA结合位点基因，仅展示5个，见表2），分析后发现gga-miR-216b可调节细胞增殖、细胞周期、凋亡和转移等过程[18]；miR-15a通过与靶基因结合降低脂肪酸氧化，间接促进鸡肌内前脂肪细胞分化[19]；miR-128-3p过表达通过调控GREM1抑制胶质瘤细胞增殖，影响细胞周期和凋亡，且LncRNA PVT1可通过miR-128-3p/GREM1轴调控BMP信号通路[20]。

表2 基因及靶向于该基因microRNATable 2 Genesand microRNAs targeting that genes

3 讨论

不同组织类型中APA亚型相对丰度存在较大差异。例如血液和睾丸倾向于表达短3´-UTR亚型，而脑组织则呈相反趋势[21]。一般情况下，增殖细胞和癌细胞中3´-UTR缩短，而在极化细胞和分化细胞中3´-UTR延长[22]。此外，不同细胞中poly（A）位点选择也存在差异。Cheng等报道合胞滋养细胞（SCT）和B细胞分化过程中存在广泛的3´-UTR缩短[23]。本研究发现前脂肪细胞来源不同，其分化过程中3´-UTR长度变化趋势也存在差异，造成差异的机制还有待进一步研究。

分析3´-UTR长度与基因本身表达水平关系发现，两者无显著相关性，因而3´-UTR缩短或延长引起microRNA结合位点缺失或增加，对基因表达量可能无影响。通过TargetScan和miRbase对部分DE-3´-UTR基因分析，发现较多microRNA靶向于该基因mRNA 3´-UTR区域。microRNA可通过与其靶基因mRNA 3´-UTR结合而影响细胞分化，且许多microRNA可直接靶作用于脂肪形成相关基因参与动物脂肪细胞增殖和分化调控，如miR-155可直接靶作用于脂肪形成的关键转录因子C/EBPβ[24]，miR-215-5p可通过靶作用于NCOA3而负调控鸡腹部脂肪形成[25]。此外，生物体内一个microRNA调控许多靶基因，且microRNA与基因mRNA靶向互作关系因组织、细胞类型和生理状态而不同[26]。因此，在前脂肪细胞分化过程中，microRNA与靶基因之间可能存在一个复杂的调控网络，下一步将重点研究脂肪形成相关靶基因鉴定和互作调控机制。

为探索前脂肪细胞不同分化阶段相似性和差异性，比较各组DEGs表达（京海黄鸡：I0vsI2、I2vsI4、I4vsI6；北京鸭：T-48vsT0、T0vsT12、T0vsT72；3T3_L1：EPvsC、Cvs10h、10hvs6d）的GO和KEGG分析结果。在京海黄鸡各组比较中发现核苷磷酸生物合成和代谢过程是I2vsI4和I4vsI6组共同的生物过程，而在I0vsI2和I4vsI6组、I0vsI2和I2vsI4组中未发现共同生物过程。表明DE-3´-UTR基因表达在分化早期（I0和I2）变化明显大于分化后期（I4和I6）。在北京鸭各组比较中发现核酸结合在T-48vsT0、T0vsT12和T0vsT72比较中存在，在T0vsT12和T0vsT72比较中有2个共同生物过程和2个共同分子功能。表明DE-3´-UTR基因表达在分化早期（T2）与分化后期（T72）差异不明显。在3T3-L1细胞各组比较中均发现细胞周期和与脂肪相关的生物过程，表明DE-3´-UTR基因在分化早期与分化后期表达具有相似性，且在分化过程中DE-3´-UTR基因表达对脂肪合成与代谢起重要调控作用。Zhang等研究与脂肪细胞分化相关的功能和通路发现，包括磷酸肌醇代谢、FoxO信号通路、TGF-beta信号通路、Wnt/β-catenin信号通路、焦点粘连、甘油磷脂新陈代谢、MAPK信号通路、脂质代谢和mTOR信号通路等[27]。GO和KEGG分析DE-3´-UTR基因发现许多功能和通路相关研究较多。例如京海黄鸡中甘油脂代谢过程、脂肪酸生物合成过程、磷酸肌醇代谢和FoxO信号通路；北京鸭中TGF-beta信号通路、焦点粘连、甘油磷脂新陈代谢和MAPK信号通路；3T3-L1细胞中MAPK活性调控、葡萄糖代谢过程、脂质代谢和mTOR信号通路。因此，这些功能和信号通路可能是进一步研究脂肪分化分子机制关键。