沙子嶺豬與大約克夏豬背最長肌miRNA-mRNA共調控網絡分析

摘 "要：微小核糖核酸（microRNA，miRNA）是一類非編碼核糖核酸（ribonuclei acid，RNA）小分子，長度約為22個核苷酸。本研究對沙子嶺豬和大約克夏豬背最長肌組織表達譜芯片和轉錄組數據進行分析，旨在揭示這兩個豬種的miRNA與信使RNA（messenger RNA，mRNA）的差異表達特征及其調控網絡。試驗通過小RNA測序共鑒定出307個miRNA，其中9個顯著差異表達（5個上調，4個下調），并分別在沙子嶺豬和大約克夏豬中鑒定出4個和5個特異性miRNA；轉錄組分析顯示鑒定出553個差異表達基因（differentially expressed genes，DEGs），其中341個上調，212個下調。整合分析發現8個差異表達miRNA（differentially expressed miRNAs，DEMs）與59個DEGs形成70對調控關系，構建了關鍵共調控網絡。基因本體論（gene ontologo，GO）富集分析表明，靶基因顯著參與肌肉收縮調控、脂肪酸氧化等生物過程；京都基因和基因組數據庫（Kyoto encyclopedia of genes and genomes，KEGG）通路分析揭示絲裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK）信號通路、FoxO（Forkhead box O，FoxO）信號通路和鐵死亡等代謝調控途徑的顯著富集，其中MAPK14基因同時參與內分泌抵抗與免疫調控，ACOX1基因在丙酸和碳代謝中發揮雙重作用。這些發現系統解析了沙子嶺豬與大約克夏豬肌肉組織在電生理特性、能量代謝和分子調控水平的關鍵差異，為解析豬種質特性形成的分子機制提供了新見解。

關鍵詞：背最長肌；miRNA芯片；Solexa測序沙子嶺豬；約克夏豬

中圖分類號：S813.1 文獻標志碼：A 文章編號：1001-0769（2025）01-0044-06

骨骼肌發育與肉質性狀的形成受多層分子調控，其中信使RNA（messenger RNA，mRNA）與微小RNA（microRNA，miRNA）的分子作用機制日益受到關注。沙子嶺豬作為中國地方豬種的典型代表，以耐粗飼、肉質細嫩著稱；而大約克夏豬作為我國引進的瘦肉型品種，則以高瘦肉率和優良生長性能而聞名于世[1]。這兩個豬種在肌纖維類型、肌內脂肪（intramuscular fat，IMF）含量及代謝通路上存在顯著差異，但其骨骼肌發育的分子調控機制尚未完全闡明。

近年來，隨著高通量測序技術的發展，發現miRNA通過靶向調控肌肉發育的相關基因，如肌球蛋白重鏈（myosin heavy chain，MYH）家族和絲裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK）家族在骨骼肌分化與脂代謝中具有關鍵作用[2]。miR-27a/b和miR-378等已被證實可通過抑制脂肪沉積影響胴體性狀[3]，而let-7 miRNA家族則參與Wnt、PI3K-Akt等信號通路的調控[4]。然而，現有研究多聚焦單一分子或組織類型，對豬種間mRNA-miRNA互作網絡的系統性分析仍較匱乏。

本研究以兩種斷奶仔豬（沙子嶺豬與大約克夏豬）背最長肌為對象，整合表達譜芯片與miRNA測序數據，通過篩選兩個品種間差異表達的mRNA與miRNA，并解析其富集的信號通路及功能關聯，構建關鍵miRNA-mRNA調控網絡，揭示品種特異性肌肉發育的分子基礎，以期為豬遺傳改良及肉質性狀分子標記的挖掘提供理論依據。

1 "材料與方法

1.1 試驗動物

試驗豬分別選自湖南正虹科技發展股份有限公司正虹原種豬場的大約克夏豬和湖南省湘潭沙子嶺豬原種場的沙子嶺豬。所選豬均為25日齡的健康去勢公豬，組內為半同胞，每個品種各3頭，飼養條件相同。屠宰后，在胸腰椎間最后一根肋骨處采集背最長肌樣本，用磷酸緩沖鹽（phosphate buffer saline，PBS）溶液沖洗后，速凍于液氮中，以備后續試驗。

1.2 總RNA提取和文庫構建

采用TRIzol法提取總RNA；通過Nanodrop 2000超微量紫外分光光度計對提取的RNA進行濃度和純度的測定；通過Agilent Bioanalyzer 2100和1.5％瓊脂糖凝膠電泳驗證RNA完整性。每個樣本取1 μL RNA進行文庫構建，測序工作由上海生物芯片有限公司采用Illumina Genome Analyzer測序平臺完成。

1.3 miRNA的鑒定和數據分析

原始數據經濾除低質量序列、ploy-N和污染序列后，保留長度為18～30 nt的序列。通過Bowtie軟件依次比對Rfam（Vol.14.1）和Repbase（Vol.24.03）數據庫，系統剔除核糖體RNA（ribosomal RNA，rRNA）、轉運RNA（transfer RNA，tRNA）、小核RNA（small nuclear RNA，snRNA）、小胞質RNA（small cytoplasmic RNA，scRNA）及重復元件等非靶序列。將過濾后的序列通過Bowtie（Vol.1.3.0）比對豬參考基因組（Sus scrofa 11.1），參數允許最大錯配數為2（不含空位錯配）。已知miRNA注釋采用miRBase（Vol.22）數據庫，新miRNA預測通過miRDeep2（Vol.2.0.1.2）完成，設置自由能閾值－20 kcal/mol且前體序列讀深≥10。各樣本miRNA表達量采用可信平臺模塊（Trusted Platform Module，TPM）算法進行標準化。基于高通量測序數據，以|log2FC（log2 fold change，log2FC）|≥1.0且P＜0.05為標準，篩選DEMs，篩選結果通過miRBase（Vol.22）注釋miRNA成熟體序列。

1.4 表達譜芯片測定與數據分析

采用Affymetrix豬全基因組表達芯片進行基因表達譜分析，每個處理組設置三個生物學重復。使用GeneChip 3' IVT Express試劑盒對6個樣本（總RNA）進行擴增、標記及純化，獲得生物素標記的互補RNA（complementary RNA，cRNA），嚴格按照試劑盒說明書操作。芯片雜交與洗滌通過GeneChip? Hybridization Wash and Stain Kit在Hybridization Oven 645雜交儀和Fluidics Station 450洗染工作站完成。采用Command Console Software 3.1控制GeneChip Scanner 3000掃描儀進行芯片掃描。原始數據經RMA算法和GeneSpringl軟件（Vol.11.0）進行背景校正及標準化處理，通過DESeq2（Vol.1.10.1）進行各樣本中mRNA的差異分析篩選，篩選標準為|log2FC|＞1.0且P＜0.05。

1.5 miRNA-mRNA調控網絡的構建

聯合TargetScan（Vol.8.0）、miRDB（Vol.6.0）數據庫對注釋后的miRNA進行靶基因預測；預測到的靶基因與轉錄組測序獲得的差異表達mRNAs（篩選標準：|log2FC|＞1.0且P＜0.05）進行維恩圖分析，獲取顯著共調控靶基因集。Cytoscape（Vol.3.7.2）構建包含DEMs及其共調控靶基因的互作網絡。

1.6 差異表達分析和富集分析

通過DESeq2（Vol.1.10.1）進行各樣本中mRNA的差異分析篩選，篩選標準為|log2FC|＞1.0且P＜0.05。參照DAVID和KOBAS 2.0數據庫對miRNA靶基因結果進行GO富集分析和KEGG通路富集分析，顯著富集的標準均為P＜0.05。

1.7 實時熒光定量逆轉錄PCR驗證分析

采用SYBR Green熒光染料法對樣本進行實時熒光定量逆轉錄PCR（real-time quantitative PCR，RT-qPCR）驗證。試驗設置3個生物學重復，以豬U6 snRNA、甲硫氨酸tRNA（Met-tRNA）及5S RNA作為內參基因，每個樣本進行3次重復。引物序列見表1和表2。采用2?ΔΔCt法計算相對表達量，組間差異通過SPSS Vol.26.0進行獨立樣本t檢驗（P＜0.05為顯著差異）。

2 "結果與分析

2.1 測序數據質量分析

原始數據（raw reads）經去接頭和低質量讀長過濾后，兩個cDNA文庫分別獲得31 862 145 " "（大約克夏豬）條和37 186 633（沙子嶺豬）條有效序列（clean reads），其中已知miRNA注釋比例分別達78.6％和77.78％（表3），表明測序數據以miRNA為主。少量序列被注釋為核仁小RNA（small nucleolar RNA，snoRNA）、核內小RNA（small nuclear RNA，snRNA）、rRNA和假基因（pseudogene）等非靶向RNA，占比較低（＜5％），表明實驗數據質量可靠。

2.2 miRNA和mRNA的差異表達分析

在沙子嶺豬和大約克夏豬的背最長肌樣本中共鑒定出307個miRNA，其中298個miRNA為大約克夏豬與沙子嶺豬共表達的，兩個品種中分別有5個（大約克夏豬）和4個（沙子嶺豬）miRNA為品種特異性表達miRNA。以大約克夏豬為對照，篩選出9個顯著差異表達的DEMs，包含5個上調的和4個下調的。

在沙子嶺豬和大約克夏豬背最長肌中的mRNA表達模式中，以大約克夏豬為對照組，共篩選出553個差異表達基因（DEGs），其中上調基因341個（62％），下調基因212個（38％）。

2.3 miRNA-mRNA調控網絡構建及其靶基因富集分析

將預測到的DEMs靶基因與轉錄組測序獲得的539個DEGs進行韋恩圖分析（Venn diagram），共獲取8個DEMs的59個基因，形成了70對調控關系（表4）。對現有的調控網絡靶基因進行GO功能富集分析，結果顯示，生物過程顯著富集于肌肉收縮調控、肌肉系統過程調節及脂肪酸氧化等通路；細胞組分主要富集于細胞器膜內在成分和離子通道復合體，其中電壓門控通道相關基因，如TF基因，呈顯著差異表達。在分子功能層面，蛋白激酶活性、電壓門控通道活性及磷酸肌醇結合等功能顯著富集，提示沙子嶺豬和大約克夏豬在肌肉電生理特性方面存在重要調控差異。KEGG通路分析發現，差異基因顯著富集于MAPK信號通路、FoxO信號通路及鐵死亡等代謝調控通路。其中，MAPK14基因同時參與內分泌抵抗和Fcε受體I信號通路調控，而ACOX1基因在丙酸代謝與碳代謝通路中均發揮重要作用。

2.4 RT-qPCR驗證

為驗證測序結果，隨機選取4個差異miRNA（ssc-miR-127、ssc-let-7e、ssc-miR-204、ssc-miR-29a）及8個DEGs（GADD45G、MYH1、MAPK14、SYMPK、SRPX、NRAP、EEA1和IL17RD）進行RT-qPCR驗證，結果顯示其表達趨勢與測序數據高度一致。

3 "討論

脂肪型豬種（中國地方品種）與瘦肉型豬種（歐洲品種）在肌肉特征上存在顯著差異[5]。以沙子嶺豬為代表的脂肪型豬種相較于瘦肉型品種具有更高的肌內脂肪[6]。肌內脂肪作為肉質的重要遺傳指標，可顯著提高肉品的嫩度、多汁性、風味和持水能力[7]，同時降低滴水損失和蒸煮損失[8]。沙子嶺豬作為華中地區代表性品種，于2014年被中國農業農村部列為重要遺傳資源保護品種，具有適應亞熱帶氣候、耐粗飼、早熟及肉質細嫩等特點[9]。

本研究發現，多個miRNA在沙子嶺豬和大約克夏豬骨骼肌中差異表達，如miR-210和miR-183可能通過調控骨骼肌肌細胞增殖與分化參與骨骼肌發育[10-11]。其中，let-7e（屬于10個亞型的let-7家族[4]）和miR-127在大約克夏豬中的表達量是沙子嶺豬的2倍，尤其是miR-127在差異miRNA中豐度最高。盡管其在豬骨骼肌生長發育中的功能尚未明確，但已有研究證實let-7家族可通過抑制血小板反應蛋白抑制劑-1和金屬蛋白酶組織抑制因子-1等抗血管生成因子促進腫瘤血管生成[12]。而與之前研究不同的是[13]，本研究發現miR-210、miR-183和miR-204在沙子嶺豬背最長肌中高表達，提示這些miRNA可能參與保育豬肌肉功能調控。

在本研究中轉錄組測序所鑒定出的539個差異基因中，共有59個差異基因受到DEMs的調控，靶基因富集分析結果表明，大約克夏與沙子嶺斷奶仔豬在肌生成和脂肪代謝方面可能存在顯著差異。GO富集分析顯示，靶基因顯著富集在涉及肌肉收縮調控、肌肉系統過程調節和脂肪酸氧化等多個肌生成相關通路。KEGG通路分析顯示，5個差異基因（ELK4/GADD45G/MKNK1/IGF1R/MAPK14）富集于MAPK信號通路。該通路不僅是骨骼肌成肌細胞分化的關鍵調控通路[14]，還可通過調控細胞連接相關通路與脂代謝通路的交互作用影響肌內脂肪沉積[15-16]。其中，Gadd45a基因編碼一種小型肌核蛋白，以刺激蛋白質分解、減少蛋白質合成、減少線粒體、抑制合成代謝信號傳導的方式改變骨骼肌基因表達，可最終導致肌纖維萎縮[17]。此外，也有研究表明，胰島素樣生長因子Ⅰ受體基因編碼的受體蛋白可結合胰島素和胰島素樣生長因子Ⅰ，通過抑制FoxO轉錄因子來維持肌肉質量和線粒體功能[18]。

4 "結論

沙子嶺豬與大約克夏豬背最長肌中miRNA及mRNA表達存在顯著差異，共篩選出9個差異miRNA（DEMs）和553個差異基因（DEGs），并構建了包含8個DEMs與59個DEGs的共調控網絡（70對互作關系）。差異表達基因進行了miRNA-mRNA調控網絡的構建以及miRNA靶基因的搜尋，對上述靶基因進行GO和KEGG富集分析，篩選出了多個與豬骨骼肌生成相關且可能受到miRNA調控的關鍵基因。其關鍵調控網絡及通路可為豬種遺傳改良及肉質性狀的分子育種提供理論依據。

參考文獻

[1] YANG H，XU X L，MA H M，et al．Integrative analysis of transcriptomics and proteomics of skeletal muscles of the Chinese indigenous Shaziling pig compared with the Yorkshire breed[J]．BMC Genetics，2016，17（1）：80.

[2] BARTEL D P．MicroRNAs genomics， biogenesis， mechanism， and function[J]．Cell，2004，116（2）：281-297.

[3] WU F，WANG F，YANG Q，et al．Upregulation of miRNA-23a-3p rescues high glucose-induced cell apoptosis and proliferation inhibition in cardiomyocytes[J]．In Vitro Cellular amp; Developmental Biology-Animal，2020，56（10）：866-877.

[4] DE ALMEIDA B C，DOS ANJOS L G，UNO M，et al． " " " " " " " " " " " " " " " " " " " " "Let-7 miRNA's expression profile and its potential prognostic role in uterine leiomyosarcoma[J]．Cells，2019，8（11）：1452.

[5] TANG Z L，LI Y，WAN P，et al．LongSAGE analysis of skeletal muscle at three prenatal stages in Tongcheng and Landrace pigs[J]．Genome Biology，2007，8（6）：R115.

[6] SUN Y M，CAI R，WANG Y Q，et al．A newly identified LncRNA LncIMF4 controls adipogenesis of porcine intramuscular preadipocyte through attenuating autophagy to inhibit lipolysis[J]．Animals，2020，10（6）：926.

[7] LIU J，NONG Q，WANG J，et al．Breed difference and regulatory role of CRTC3 in porcine intramuscular adipocyte[J]．Animal Genetics，2020，51（4）：521-530.

[8] GERBENS F，VERBURG F J，VAN MOERKERK H T，et al．Associations of heart and adipocyte fatty acid-binding protein gene expression with intramuscular fat content in pigs[J]．Journal of Animal Science，2001，79（2）：347-354.

[9] XING X W，XUE L Q，HUANG S Q，et al．Research of porcine endogenous retrovirus in Shaziling pigs[J]．Yi Chuan，2006，28（7）：799-804.

[10] ISMAEEL A，FLETCHER E，MISERLIS D，et al．Skeletal muscle MiR-210 expression is associated with mitochondrial function in peripheral artery disease patients[J]．Translational Research，2022，246：66-77.

[11] ZHANG L Y，ZHONG D D，YAO C X，et al．Buffalo bbu-miR-493-5p promotes myoblast proliferation and differentiation[J]．Animals，2024，14（4）：533.

[12] OTSUKA M，ZHENG M，HAYASHI M，et al．Impaired microRNA processing causes corpus luteum insufficiency and infertility in mice[J]．The Journal of Clinical Investigation，2008，118（5）：1944-1954.

[13] CHEN C，DENG B，QIAO M，et al．Solexa sequencing identification of conserved and novel microRNAs in backfat of Large White and Chinese Meishan pigs[J]．PLoS One，2012，7（2）：e31426.

[14] KEREN A，TAMIR Y，BENGAL E．The p38 MAPK signaling pathway： A major regulator of skeletal muscle development[J]．Molecular and Cellular Endocrinology，2006，252（1/2）：224-230.

[15] KOKTA T A，DODSON M V，GERTLER A，et al．Intercellular signaling between adipose tissue and muscle tissue[J]．Domestic Animal Endocrinology，2004，27（4）：303-331.

[16] CUI H X，LIU R R，ZHAO G P，et al．Identification of differentially expressed genes and pathways for intramuscular fat deposition in pectoralis major tissues of fast-and slow-growing chickens[J]．BMC Genomics，2012，13：213.

[17] BONGERS K S，FOX D K，EBERT S M，et al．Skeletal muscle denervation causes skeletal muscle atrophy through a pathway that involves both Gadd45a and HDAC4[J]．American Journal of Physiology．Endocrinology and Metabolism，2013，305（7）：E907- E915.

[18] BHARDWAJ G，PENNIMAN C M，KLAUS K，et al．Transcriptomic regulation of muscle mitochondria and calcium signaling by insulin/IGF-1 receptors depends on FoxO transcription factors[J]．Frontiers in Physiology，2022，12：779121.