收稿日期：2023-09-18

基金項(xiàng)目：國(guó)家兔產(chǎn)業(yè)技術(shù)體系項(xiàng)目（CARS-43-D-1）；四川省財(cái)政運(yùn)行專項(xiàng)（SASA2024CZYX004）；四川省“十四五”育種攻關(guān)項(xiàng)目（2021YFYZ0033）；四川省科研院所科技成果轉(zhuǎn)化項(xiàng)目（2023N22J0001）

作者簡(jiǎn)介：曾建紅（1998-），女，四川綿陽(yáng)人，碩士，主要從事家兔繁殖與飼養(yǎng)管理。（E-mail）1762725547@qq.com

通訊作者：鄺良德，（E-mail）215640832@qq.com

摘要： 為明確蜀興1號(hào)肉兔（SX）與伊拉兔（IRA）肉質(zhì)差異及差異形成分子機(jī)制，本研究通過(guò)采集2個(gè)品種兔背最長(zhǎng)肌樣品，進(jìn)行肉質(zhì)性狀測(cè)定及轉(zhuǎn)錄組測(cè)序和代謝組分析，篩選差異表達(dá)基因和差異代謝物并進(jìn)行功能富集分析。結(jié)果表明，蜀興1號(hào)肉兔的肌纖維直徑和滴水損失率顯著低于伊拉兔，而肌纖維密度、熟肉率和肌內(nèi)脂肪含量顯著高于伊拉兔。蜀興1號(hào)肉兔與伊拉兔轉(zhuǎn)錄組中共篩選出81個(gè)差異表達(dá)基因，其中蜀興1號(hào)肉兔中51個(gè)基因上調(diào)，30個(gè)基因下調(diào)。SMTNL1、PM20D2和EDN1等可能是導(dǎo)致兩種兔肉品質(zhì)差異的基因。差異表達(dá)基因顯著富集于cAMP信號(hào)傳導(dǎo)途徑。揮發(fā)性代謝組學(xué)比較共得到12種差異顯著代謝物，均在蜀興1號(hào)肉兔中上調(diào)，其中，2-十一烯醛、4-乙基辛酸、（E）-2-壬烯醛、鳥氨酸和十一醛等代謝物與肉質(zhì)風(fēng)味正相關(guān)，揮發(fā)性代謝物顯著富集于2-氧代羧酸代謝、ABC轉(zhuǎn)運(yùn)蛋白和精氨酸生物合成等KEGG代謝通路。廣泛靶向代謝組學(xué)比較分析共得到15種差異顯著代謝物，其中蜀興1號(hào)肉兔中7種代謝物上調(diào)，8種下調(diào)，2個(gè)品種兔肉中γ-L-谷氨酸-L-谷氨酰胺、L-谷氨酸-L-谷氨酰胺和溶血磷脂酰膽堿（LPC）等肉質(zhì)風(fēng)味相關(guān)物質(zhì)的含量差異較大，廣泛靶向代謝物顯著富集于咖啡因代謝和晝夜節(jié)律夾帶等KEGG代謝通路。本研究得到的差異表達(dá)基因和差異代謝物可為進(jìn)一步的優(yōu)質(zhì)兔養(yǎng)殖和優(yōu)質(zhì)兔選育提供參考依據(jù)。

關(guān)鍵詞： 轉(zhuǎn)錄組學(xué)；代謝組學(xué)；蜀興1號(hào)肉兔；伊拉兔；肉質(zhì)

中圖分類號(hào)： S829.1"" 文獻(xiàn)標(biāo)識(shí)碼： A"" 文章編號(hào)： 000-4440（2024）09-1689-12

Molecular mechanism of transcriptomic and metabolomic differences between the longissimus dorsi muscle of Shuxing No.1 rabbit and Ira rabbit

ZENG Jianhong ZHENG Yucai LI Congyan GUO Zhiqiang YANG Rui ZHENG Jie LI Yuying REN Yongjun LEI Min XIE Xiaohong KUANG Liangde ，3

（1.Sichuan Animal Science Academy， Chengdu 610066， China;2.College of Animal and Veterinary Sciences， Southwest Minzu University， Chengdu 610041， China;3.Animal Genetics and Breeding Key Laboratory of Sichuan Province， Chengdu 610066， China）

Abstract： In order to clarify the difference of meat quality between Shuxing No.1 rabbit （SX） and Ira rabbit （IRA） and the molecular mechanism of the difference formation， the longissimus dorsi muscle samples of two varieties of rabbits were collected for meat quality trait determination， transcriptome sequencing and metabolomics analysis， and differentially expressed genes and differentially expressed metabolites were screened and their functional enrichment analysis was performed. The results showed that the muscle fiber diameter and drip loss rate of SX were significantly lower than those of IRA， while the muscle fiber density， cooked meat rate and intramuscular fat content were significantly higher than those of IRA. A total of 81 differentially expressed genes were screened in the transcriptome of SX and IRA. Among them， 51 genes were up-regulated and 30 genes were down-regulated in SX. SMTNL1， PM20D2 and EDN1 might be the genes that lead to the difference in meat quality between SX and IRA. The differentially expressed genes were significantly enriched in the cAMP signaling pathway. A total of 2 significantly different metabolites were obtained by volatile metabolomics comparison， all of which were up-regulated in SX. Among them， metabolites such as 2-undecenal， 4-ethyloctanoic acid， （E）-2-nonenal， ornithine and undecanal were positively correlated with meat flavor. The volatile metabolites were significantly enriched in KEGG metabolic pathways such as 2-oxocarboxylic acid metabolism， ABC transporter and arginine biosynthesis. A total of 5 significantly different metabolites were obtained by extensive targeted metabolomics comparative analysis. Among them， seven metabolites were up-regulated and eight metabolites were down-regulated in SX. The contents of meat flavor-related substances such as γ-L-glutamic acid-L-glutamine， L-glutamic acid-L-glutamine and lysophosphatidylcholine （LPC） in SX and IRA were significantly different. Widely targeted metabolites were significantly enriched in KEGG metabolic pathways such as caffeine metabolism and circadian rhythm entrainment. The differentially expressed genes and differential metabolites obtained in this study can provide a reference for further high-quality rabbit breeding.

Key words： transcriptomics;metabolomics;Shuxing No.1 rabbit;Ira rabbit;meat quality

中國(guó)是世界上第一肉兔養(yǎng)殖大國(guó)，2022年兔肉產(chǎn)量5.18×10 5t，兔業(yè)產(chǎn)值超過(guò)3.00×10 0元。優(yōu)質(zhì)畜禽資源的培育與養(yǎng)殖是畜禽產(chǎn)業(yè)健康發(fā)展的基礎(chǔ)，同時(shí)符合消費(fèi)者的需要。目前中國(guó)兔產(chǎn)業(yè)總體呈上升趨勢(shì)，但由于對(duì)種兔重要性認(rèn)知不夠，導(dǎo)致中國(guó)優(yōu)良種兔的繁育體系不健全、優(yōu)質(zhì)兔資源缺乏及品種退化，因此長(zhǎng)期以來(lái)只能通過(guò)從國(guó)外引進(jìn)優(yōu)良品種以維持兔產(chǎn)業(yè)的發(fā)展，這不但增加了兔養(yǎng)殖成本，還制約著兔產(chǎn)業(yè)的發(fā)展[1]。伊拉兔是由法國(guó)歐洲兔業(yè)公司在20世紀(jì)70年代末培育成的肉兔配套系，具有早期生長(zhǎng)速度快、抗病力強(qiáng)、出肉率高等特點(diǎn)，因其適應(yīng)中國(guó)大部分地區(qū)的氣候和飼養(yǎng)條件，被廣泛飼養(yǎng)和銷售[2]。近年來(lái)，中國(guó)的優(yōu)質(zhì)兔育種取得了一定的成績(jī)，培育出了一些優(yōu)質(zhì)兔品種。蜀興1號(hào)肉兔是四川省畜牧科學(xué)研究院針對(duì)中國(guó)西南地區(qū)特殊的兔肉消費(fèi)需求和粗放的飼養(yǎng)管理?xiàng)l件，利用齊興肉兔、歐洲大白兔和齊卡新西蘭白兔3個(gè)品種（系）肉兔培育而成的優(yōu)質(zhì)肉兔配套系，具有繁殖性能好、適應(yīng)性強(qiáng)、上市期早、屠宰率高、耗料少等特點(diǎn)，于2020年12月通過(guò)國(guó)家畜禽遺傳資源委員會(huì)審定[3]。

基于轉(zhuǎn)錄組測(cè)序技術(shù)和代謝組學(xué)技術(shù)分析畜禽產(chǎn)品品質(zhì)差異的研究已得到廣泛應(yīng)用。牛、羊、鴨、雞等畜禽產(chǎn)品中品質(zhì)相關(guān)基因的挖掘、功能注釋及代謝通路已有大量研究[4-7]，雞肉、鴨肉、豬肉等畜禽產(chǎn)品中與肉質(zhì)風(fēng)味相關(guān)的代謝物及前體物質(zhì)篩選亦有初步研究[8-11]。截止目前，基于轉(zhuǎn)錄組測(cè)序和代謝組學(xué)技術(shù)進(jìn)行兔肉品質(zhì)的比較及品質(zhì)差異形成的分子機(jī)制研究還不多見。本研究以蜀興1號(hào)肉兔和伊拉兔為研究對(duì)象，利用轉(zhuǎn)錄組測(cè)序和代謝組學(xué)技術(shù)，分析蜀興1號(hào)肉兔和伊拉兔差異表達(dá)基因及其代謝通路，為優(yōu)質(zhì)兔的養(yǎng)殖和分子育種提供參考。

1 材料與方法

1.1 試驗(yàn)材料及肉質(zhì)性能測(cè)定

試驗(yàn)動(dòng)物為28日齡商品代伊拉兔（IRA）（購(gòu)買于重慶阿興記原種兔養(yǎng)殖基地）和28日齡蜀興1號(hào)肉兔（SX）（四川省畜牧科學(xué)研究院繁育），于四川省畜牧科學(xué)研究院肉兔科研基地相同條件飼養(yǎng)42 d，選取蜀興1號(hào)肉兔和伊拉兔各30只（公兔、母兔各15只）進(jìn)行屠宰，將兔背最長(zhǎng)肌組織按照檢測(cè)要求進(jìn)行分割，冷藏保存。屠宰后45 min和24 h，利用Testo 205 pH測(cè)量?jī)x（德國(guó)Testo公司產(chǎn)品）測(cè)定肌肉的pH值，利用美能達(dá)色差儀CR-10（日本美能達(dá)公司產(chǎn)品）測(cè)定兔肉的亮度值（L *）、紅綠度（a *）、黃藍(lán)度（b *）。參照劉浪等[12]的方法測(cè)定兔背最長(zhǎng)肌的肌纖維直徑和密度、熟肉率、滴水損失率、脂肪含量、蛋白質(zhì)含量、水分含量和灰分含量。使用IBM SPSS統(tǒng)計(jì)軟件進(jìn)行分析差異顯著性。試驗(yàn)兔屠宰30 min內(nèi)，將兔背最長(zhǎng)肌組織按照測(cè)序要求分割放置于凍存管中，液氮速凍后置于-80 ℃超低溫冰箱保存。隨機(jī)選取蜀興1號(hào)肉兔和伊拉兔各8只（公、母各4只）由武漢邁特維爾生物科技有限公司完成轉(zhuǎn)錄組學(xué)和代謝組學(xué)測(cè)序。

1.2 主要儀器與試劑

試驗(yàn)用儀器主要有美國(guó)Thermo Fisher Scientific公司生產(chǎn)的Qubit 2.0熒光分光光度計(jì)、Thermo Sorvall ST 8R高速冷凍離心機(jī)和核酸蛋白檢測(cè)儀，美國(guó)Bio-Rad公司生產(chǎn)的CFX96 Touch Real-Time PCR Detection System、POWER/PAC 3000電泳儀及VersaDocTM凝膠成像系統(tǒng)，美國(guó)Sciex公司生產(chǎn)的ExionLC AD超高效液相色譜系統(tǒng)，美國(guó)Agilent公司生產(chǎn)的Agilent 2100 生物分析儀、7890B-7000D氣相色譜-質(zhì)譜聯(lián)用儀、毛細(xì)管色譜柱DB-5MS（規(guī)格30 m×0.25 mm×0.25 μm），美國(guó)Waters公司生產(chǎn)的廣泛靶向色譜柱Waters ACQUITY UPLC HSS T3 C18（規(guī)格1.8 μm，2.1 mm×100.0 mm）；德國(guó)Implen公司生產(chǎn)的Nano Photometer 超微量分光光度計(jì)；德國(guó)Retsch公司生產(chǎn)的MM400球磨儀，瑞士CTC Analytics AG公司生產(chǎn)的SPME Arrow固相微萃取裝置、Fiber Conditioning Station老化裝置、Agitator樣品加熱箱。

主要試劑有美國(guó)Illumina公司生產(chǎn)的NEBNextUltraTM RNA Library Prep Kit試劑盒，美國(guó)Thermo Fisher Scientific公司生產(chǎn)的Thermo K1622逆轉(zhuǎn)錄試劑盒、TRIzolReagent，日本TaKaRa公司生產(chǎn)的TB Green Premix Ex Taq Ⅱ、DL2000 DNA Marker，北京金沙生物技術(shù)有限公司生產(chǎn)的GS-GelRed核酸凝膠染料，美國(guó)Merck公司生產(chǎn)的正己烷、甲醇、乙腈，上海Aladdin公司生產(chǎn)的甲酸、甲酸銨、氨水，廣州賽國(guó)生物科技公司生產(chǎn)的氯仿、異丙醇、無(wú)水乙醇、氯化鈉（分析純）。

1.3 兔肉總RNA提取及轉(zhuǎn)錄組測(cè)序

利用TRIzol法提取背最長(zhǎng)肌組織中的總RNA，用核酸蛋白檢測(cè)儀檢測(cè)RNA的純度，用Qubit 2.0熒光分光光度計(jì)、Agilent 2100 生物分析儀檢測(cè)RNA濃度和完整性，使用NEBNext UltraTM RNA Library Prep Kit試劑盒，通過(guò)磁珠富集mRNA，并合成cDNA，再使用Qubit2.0熒光分光光度計(jì)和Agilent 2100生物分析儀對(duì)cDNA的濃度和完整性進(jìn)行檢測(cè)，最后用AMPure XP磁珠篩選出200 bp左右的cDNA，進(jìn)行PCR富集得到最終的cDNA文庫(kù)。文庫(kù)構(gòu)建完成后，通過(guò)qPCR方法對(duì)文庫(kù)有效濃度進(jìn)行準(zhǔn)確定量，對(duì)文庫(kù)質(zhì)量進(jìn)行檢驗(yàn)，庫(kù)檢合格后，通過(guò)Illumina測(cè)序平臺(tái)進(jìn)行高通量測(cè)序。使用Fastq軟件對(duì)原始數(shù)據(jù)進(jìn)行質(zhì)量檢驗(yàn)，去除A、T堿基以及G、C堿基分離的低質(zhì)量數(shù)據(jù)，得到高質(zhì)量數(shù)據(jù)，并對(duì)高質(zhì)量數(shù)據(jù)中的G、C堿基含量進(jìn)行檢測(cè)。使用String Tie軟件進(jìn)行新基因預(yù)測(cè)，使用Feature Count軟件對(duì)每個(gè)樣本的基因數(shù)據(jù)進(jìn)行統(tǒng)計(jì)并合并所有樣品的基因計(jì)數(shù)結(jié)果[13-14]。采用FPKM（Fragments Per Kilobase of exon model per Million mapped fragments）值衡量基因表達(dá)水平，反映轉(zhuǎn)錄本表達(dá)水平。對(duì)于無(wú)生物學(xué)重復(fù)的樣品使用edgeR差異分析軟件獲得組間差異表達(dá)基因集[15]。利用HISAT2軟件將高質(zhì)量數(shù)據(jù)與參考基因組Oryctolagus_cuniculus.OryCun2.0.dna.toplevel.fa.gz（http：//ftp.ensembl.org/pub/release-105/gff3/oryctolagus_cuniculus/）進(jìn)行序列比對(duì)后根據(jù)基因的表達(dá)量進(jìn)行差異表達(dá)分析，當(dāng)|log2FC|≥1（FC為差異倍數(shù)），且FDR（False discovery rate，錯(cuò)誤發(fā)現(xiàn)率）lt;0.05時(shí)，認(rèn)為該基因在2個(gè)兔品種間呈差異表達(dá)[16]。

1.4 差異表達(dá)基因的qPCR驗(yàn)證分析

為驗(yàn)證RNA-seq所得到的差異表達(dá)基因的可靠性，選取表達(dá)差異較大的4個(gè)基因（MCAM、SLC7A6、LCMT2和CALB2），以家兔GAPDH和18SrRNA作為內(nèi)參基因，利用qPCR技術(shù)對(duì)4個(gè)基因的表達(dá)模式進(jìn)行驗(yàn)證?；蛞锶绫?所示。

1.5 代謝組學(xué)測(cè)序

1.5.1 GC-MS揮發(fā)性代謝組學(xué)測(cè)序 在頂空瓶中添加0.2 g兔背最長(zhǎng)肌樣品、0.2 g NaCl粉末和10 μL內(nèi)標(biāo)溶液，先在60 ℃下靜置5 min，然后插入萃取頭并在250 ℃條件下解析5 min，使用SPME Arrow全自動(dòng)頂空固相微萃取裝置進(jìn)行樣本萃取。色譜條件：采用DB-5MS毛細(xì)管柱（30 m×0.25 mm×0.25 μm），載氣為高純氦氣，流速為1.2 mL/min，進(jìn)樣口溫度設(shè)定為250 ℃，溶劑延遲3.5 min。升溫程序設(shè)置為開始時(shí)保持40 ℃ 3.5 min，以10 ℃/min的速率升至100 ℃，再以7 ℃/min的速率升至180 ℃，最后以25 ℃/min的速率升至280 ℃并保持5 min。質(zhì)譜條件為：采用電子轟擊離子源，溫度為230 ℃，四級(jí)桿溫度150 ℃，質(zhì)譜接口溫度280 ℃，電子能量設(shè)定為70 eV，掃描方式為離子檢測(cè)模式（SIM），定性定量離子精準(zhǔn)掃描。每種化合物分別選擇2～3個(gè)定性離子和1個(gè)定量離子。然后按照質(zhì)譜出峰順序分時(shí)段檢測(cè)所有離子，并選取其中的定量離子進(jìn)行積分和校正工作，得到各種化合物的含量。

1.5.2 廣泛靶向代謝組學(xué) 將20 mg樣品放到離心管中，并加入5 mm鋼珠，先用MM400球磨儀以30 Hz的頻率均勻攪拌20 s，然后4 ℃條件下以3 000 r/min的速度離心30 s，離心后加入70%甲醇400 μL，振蕩5 min后放于冰上靜置15 min，再在4 ℃條件下以12 000 r/min離心10 min，將200 μL上清液移到進(jìn)樣瓶中，用于上機(jī)分析。色譜條件為：流動(dòng)相A為超純水，流動(dòng)相B為乙腈；洗脫梯度：0～11.9 min水與乙腈體積比為95∶5，12.0～12.1 min水與乙腈體積比為10∶90，12.2～14.0 min水與乙腈體積比為95∶5；流速0.4 mL/min；柱溫40 ℃；進(jìn)樣量為2 μL。質(zhì)譜條件為：電噴霧離子源溫度500 ℃，正電壓5 500 V、負(fù)電壓-4 500 V；離子源氣體設(shè)定為氣體Ⅰ379.3 kPa、氣體Ⅱ413.8 kPa、氣簾氣172.4 kPa。基于自建靶向標(biāo)準(zhǔn)品數(shù)據(jù)庫(kù)MWDB，通過(guò)檢測(cè)物質(zhì)的子、母離子和保留時(shí)間對(duì)二級(jí)譜數(shù)據(jù)進(jìn)行準(zhǔn)確判斷。利用MultiQuant軟件對(duì)不同樣本中相同代謝物的質(zhì)譜峰值進(jìn)行積分校正，并利用多反應(yīng)監(jiān)測(cè)模式進(jìn)行代謝物的定量分析。

1.6 轉(zhuǎn)錄組學(xué)和代謝組學(xué)的數(shù)據(jù)分析

將所有比較組的差異基因合并后取并集作為差異基因集進(jìn)行層次聚類分析，使用Z分?jǐn)?shù)對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理并繪制聚類熱圖。采用基因功能國(guó)際標(biāo)準(zhǔn)分類體系GO和綜合性數(shù)據(jù)庫(kù)KEGG對(duì)差異基因進(jìn)行分析，進(jìn)一步解讀基因的功能。

利用Mass Hunter軟件處理GC-MS揮發(fā)性代謝組學(xué)獲得的代謝物數(shù)據(jù)，利用Analyst .6.3軟件處理廣泛靶向代謝組學(xué)獲得的質(zhì)譜數(shù)據(jù)。利用主成分分析（PCA）和正交偏最小二乘法-判別分析（OPLS-DA）評(píng)估2個(gè)兔品種間的代謝物差異。基于變量重要性投影值（VIP）≥1.0和|log2FC|≥1，篩選出蜀興1號(hào)肉兔和伊拉兔背最長(zhǎng)肌的差異代謝物，利用綜合性數(shù)據(jù)庫(kù)KEGG進(jìn)行差異代謝物的代謝途徑富集分析。

2 結(jié)果與分析

2.1 肉質(zhì)性狀

蜀興1號(hào)肉兔和伊拉兔屠宰后背最長(zhǎng)肌組織的肉質(zhì)性狀如表2所示。從表中可知，2個(gè)品種兔屠宰后45 min和24 h的pH值和肉色指標(biāo)（亮度值L *、紅度值a *和黃度值b *）無(wú)顯著差異；隨著貯藏時(shí)間的增加，兔肉的pH值下降，酸性增強(qiáng)，而L *、a *和b *均呈增加趨勢(shì)。蜀興1號(hào)肉兔的肌纖維直徑和滴水損失率顯著低于伊拉兔，而肌纖維密度和熟肉率顯著高于伊拉兔。兔肉常規(guī)組分中，蜀興1號(hào)肉兔的肌內(nèi)脂肪含量顯著高于伊拉兔，而2個(gè)品種兔肉的含水率、蛋白質(zhì)含量、灰分含量差異不顯著。

2.2 轉(zhuǎn)錄組學(xué)測(cè)序統(tǒng)計(jì)及質(zhì)量評(píng)估

本研究建立的16個(gè)測(cè)序文庫(kù)共獲得約8.5×10 8 bp的原始數(shù)據(jù)，質(zhì)量控制后，約97%的原始數(shù)據(jù)達(dá)到高質(zhì)量數(shù)據(jù)標(biāo)準(zhǔn)。Q20（質(zhì)量值≥20的堿基所占百分比）均在96%以上，Q30（質(zhì)量值≥30的堿基所占百分比）均在90%以上，說(shuō)明轉(zhuǎn)錄組學(xué)測(cè)序結(jié)果質(zhì)量較好。G+C堿基含量在50.02%～56.59%，符合理論分布比例，說(shuō)明本研究測(cè)序結(jié)果是可靠的。16個(gè)測(cè)序文庫(kù)質(zhì)量控制后的數(shù)據(jù)與參考基因組的比對(duì)率為79.95%～87.49%（表3）。

2.3 蜀興1號(hào)肉兔和伊拉兔差異表達(dá)基因的篩選

蜀興1號(hào)肉兔和伊拉兔背最長(zhǎng)肌共有81個(gè)差異表達(dá)基因。其中，蜀興1號(hào)肉兔中表達(dá)上調(diào)的差異基因51個(gè)，表達(dá)下調(diào)的差異基因30個(gè)（圖1）。蜀興1號(hào)肉兔中差異倍數(shù)較大的表達(dá)上調(diào)基因包括SLC7A6、MCAM、LARS2和EDN1等，表達(dá)下調(diào)基因包括CALB2、LCMT2、SMTNL1和PM20D2等。

蜀興1號(hào)肉兔和伊拉兔背最長(zhǎng)肌差異表達(dá)基因?qū)哟尉垲惤Y(jié)果如圖2所示。從圖中可以看出，16個(gè)轉(zhuǎn)錄本按兔品種聚為2類，說(shuō)明這些聚類基因可能具有相似的功能或者處于同一代謝通路。

2.4 蜀興1號(hào)肉兔和伊拉兔差異表達(dá)基因的功能富集分析

蜀興1號(hào)肉兔和伊拉兔81個(gè)差異表達(dá)基因共富集到758條GO條目中，包括79個(gè)細(xì)胞組分條目、578個(gè)生物進(jìn)程條目、101個(gè)分子功能條目。主要功能條目中富集到的基因數(shù)如圖3所示。細(xì)胞過(guò)程、代謝過(guò)程、生物過(guò)程調(diào)控和生物調(diào)控等生物進(jìn)程中聚集的基因數(shù)較多；分子功能中以結(jié)合和催化活性等功能富集的基因數(shù)量較多；細(xì)胞組分中以細(xì)胞解剖實(shí)體、蛋白質(zhì)復(fù)合體等功能富集的基因數(shù)較多。81個(gè)差異表達(dá)基因KEGG富集分析后共得到61條KEGG通路，其中富集最顯著的20個(gè)通路如圖4所示。從圖4可以看出，富集到cAMP信號(hào)傳導(dǎo)途徑的差異表達(dá)基因數(shù)量最豐富、顯著度亦高，富集到類固醇激素生物合成和RNA降解途徑的差異表達(dá)基因數(shù)量和富集顯著度略低。

2.5 蜀興1號(hào)肉兔和伊拉兔差異表達(dá)基因的qPCR驗(yàn)證

基于qPCR和RNA-seq得到的MCAM、SLC7A6、LCMT2和CALB2等4個(gè)差異表達(dá)基因表達(dá)特征如圖5所示。從圖中可以看出，qPCR和RNA-seq得到的MCAM、SLC7A6、LCMT2和CALB2等4個(gè)基因的表達(dá)趨勢(shì)基本一致，說(shuō)明RNA-seq測(cè)序結(jié)果準(zhǔn)確可靠，能夠體現(xiàn)蜀興1號(hào)肉兔和伊拉兔背最長(zhǎng)肌中差異表達(dá)基因的表達(dá)特征。

2.6 代謝產(chǎn)物分析與鑒定

基于GC-MS揮發(fā)性代謝組學(xué)分析方法，蜀興1號(hào)肉兔和伊拉兔中共檢測(cè)到368種代謝物，其中12種代謝物存在顯著差異，主要包括2-十一烯醛、5-氨基-2，2，4-三甲基環(huán)戊烷甲胺和4-乙基辛酸等，且都在蜀興1號(hào)肉兔中顯著上調(diào)（表4）?；趶V泛靶向代謝組學(xué)分析方法，在蜀興1號(hào)肉兔和伊拉兔中共檢測(cè)到734種代謝物，15種代謝物存在顯著差異。其中，γ-L-谷氨酸-L-谷氨酰胺、L-谷氨酸-L-谷氨酰胺、褪黑激素等7種代謝物在蜀興1號(hào)肉兔中顯著上調(diào)，溶血磷脂酰膽堿、DL-泛酰醇、4-羥基苯甲酸丙酯等8種代謝物在蜀興1號(hào)肉兔中顯著下調(diào)。蜀興1號(hào)肉兔中γ-L-谷氨酸-L-谷氨酰胺含量和L-谷氨酸-L-谷氨酰胺含量約為伊拉兔的14.90倍、褪黑激素含量約為伊拉兔的9.34倍。伊拉兔中溶血磷脂酰膽堿含量約為蜀興1號(hào)肉兔的14.98倍（表5）。γ-L-谷氨酸-L-谷氨酰胺、L-谷氨酸-L-谷氨酰胺、褪黑激素和溶血磷脂酰膽堿可能是造成蜀興1號(hào)肉兔和伊拉兔品質(zhì)差異的關(guān)鍵代謝物。

2.7 差異代謝物KEGG富集分析

根據(jù)KEGG數(shù)據(jù)庫(kù)，基于GC-MS揮發(fā)性代謝組學(xué)技術(shù)得到的蜀興1號(hào)肉兔和伊拉兔368種代謝物共富集到8種代謝通路，包括2-氧代羧酸代謝、ABC轉(zhuǎn)運(yùn)蛋白、精氨酸生物合成、氨基酸的生物合成、D-氨基酸代謝、谷胱甘肽代謝、精氨酸和脯氨酸代謝等（圖6A）。基于廣泛靶向代謝組學(xué)技術(shù)得到的734種代謝物共富集到咖啡因代謝、晝夜節(jié)律夾帶、色氨酸代謝、泛酸和輔酶A生物合成、神經(jīng)活性配體-受體相互作用等代謝通路，其中，最可能影響蜀興1號(hào)肉兔和伊拉兔背最長(zhǎng)肌代謝物變化的代謝途徑為咖啡因代謝（圖6B）。

3 討論與結(jié)論

肉品質(zhì)和風(fēng)味是評(píng)價(jià)畜禽肉質(zhì)優(yōu)劣的重要指標(biāo)，同時(shí)也是畜禽品種改良的關(guān)鍵指標(biāo)之一[17]。肌內(nèi)脂肪、氨基酸、核苷酸和有機(jī)酸等物質(zhì)與肉的風(fēng)味密切相關(guān)，含量的不同會(huì)引起肉品質(zhì)和風(fēng)味產(chǎn)生差異[18]。品種、年齡和營(yíng)養(yǎng)狀況等多種因素都會(huì)導(dǎo)致畜禽產(chǎn)品肌纖維特性的差異。通常情況下，肌肉中肌纖維密度越高、肌纖維直徑越小，其細(xì)嫩度越佳[19]。本研究結(jié)果表明，蜀興1號(hào)肉兔和伊拉兔的肉品質(zhì)指標(biāo)存在一定程度的差異。蜀興1號(hào)肉兔的肌纖維密度顯著高于伊拉兔，肌纖維直徑顯著低于伊拉兔。這說(shuō)明，蜀興1號(hào)肉兔的肉質(zhì)比伊拉兔更加細(xì)嫩。

轉(zhuǎn)錄組學(xué)研究可揭示畜禽生長(zhǎng)和品質(zhì)差異的關(guān)鍵基因[20-26]，對(duì)優(yōu)質(zhì)畜禽資源的培育具有較好的指導(dǎo)意義。本研究通過(guò)轉(zhuǎn)錄組測(cè)序技術(shù)，共篩選出蜀興1號(hào)肉兔和伊拉兔背最長(zhǎng)肌中81個(gè)差異表達(dá)基因，其中，51個(gè)基因在蜀興1號(hào)肉兔中表達(dá)上調(diào)，30個(gè)下調(diào)。SMTNL1、PM20D2和EDN1等基因在2個(gè)品種兔之間表達(dá)差異顯著，而這些基因已被證實(shí)能調(diào)控畜禽產(chǎn)品品質(zhì)。SMTNL1能調(diào)節(jié)兔平滑肌收縮力，加強(qiáng)平滑肌和骨骼肌對(duì)高血壓、妊娠和運(yùn)動(dòng)訓(xùn)練的適應(yīng)，還可能參與調(diào)控細(xì)胞骨架的重組和細(xì)胞收縮[27]。SMTNL1和孕酮受體之間的相互作用能改變收縮蛋白和代謝蛋白的表達(dá)，促進(jìn)小鼠懷孕期間骨骼肌纖維轉(zhuǎn)換，SMTNL1的缺失會(huì)降低小鼠的代謝效率和葡萄糖耐受量[28]。SMTNL1蛋白還能夠與肌球蛋白磷酸酶等收縮調(diào)節(jié)因子相互作用，在肌肉組織和類固醇激素敏感組織中表達(dá)[29-30]。這些研究結(jié)果都說(shuō)明SMTNL1在肌肉生長(zhǎng)和能量代謝中具有重要作用。PM20D2是1種含有M20酶結(jié)構(gòu)域的蛋白質(zhì)，而M20酶結(jié)構(gòu)域蛋白質(zhì)家族主要與脂質(zhì)代謝相關(guān)，在能量代謝和脂質(zhì)代謝中起重要作用[31]。PM20D2還能阻止β-丙氨酰賴氨酸等異常二肽的積累，從而有利于肌肽和高肌肽的合成[32]。本研究中差異表達(dá)基因的KEGG富集分析結(jié)果表明，富集到cAMP信號(hào)傳導(dǎo)途徑的差異表達(dá)基因數(shù)量最多，且顯著度最高，而cAMP信號(hào)傳導(dǎo)途徑通路已被證實(shí)是與肉品質(zhì)調(diào)控和風(fēng)味物質(zhì)代謝有關(guān)的通路，富集到cAMP信號(hào)傳導(dǎo)途徑通路的EDN1可通過(guò)自分泌和旁分泌作用促進(jìn)血管平滑肌細(xì)胞的增殖或肥大，從而影響平滑肌細(xì)胞的收縮[33]。因此，本研究篩選出的SMTNL1、PM20D2和EDN1差異表達(dá)基因可能是導(dǎo)致蜀興1號(hào)肉兔和伊拉兔肉質(zhì)性能差異的關(guān)鍵基因。

脂肪酸、氨基酸和醛類等代謝物是影響動(dòng)物肉質(zhì)和風(fēng)味的重要因素[34]。十一烯醛是牛脂酶解-輕度熱氧化過(guò)程中產(chǎn)生的獨(dú)特產(chǎn)物，是牛肉加熱過(guò)程中產(chǎn)生的香氣物質(zhì)[35]，（E）-2-壬烯醛、2-十一烯醛等代謝物是灘羊肉的關(guān)鍵香氣化合物[36]，4-乙基辛酸則是羊脂中的主要風(fēng)味活性物質(zhì)[37]。本研究通過(guò)GC-MS揮發(fā)性代謝組學(xué)測(cè)序技術(shù)，共篩選出12個(gè)差異顯著的代謝物，包括2-十一烯醛、4-乙基辛酸、（E）-2-壬烯醛和十一醛等，這些物質(zhì)均在蜀興1號(hào)肉兔肉中顯著上調(diào)。通過(guò)廣泛靶向代謝組學(xué)測(cè)序技術(shù)，共篩選出差異顯著的代謝物15個(gè)，其中蜀興1號(hào)肉兔肉中上調(diào)代謝物7個(gè)，下調(diào)代謝物8個(gè)。差異代謝物大多顯著富集在咖啡因代謝通路。在篩選出的差異代謝物中，肉堿能將長(zhǎng)鏈脂肪酸從細(xì)胞質(zhì)輸送到線粒體基質(zhì)中，提高脂肪酸的轉(zhuǎn)化效率[38]；咖啡因可以防止脂肪在細(xì)胞中過(guò)度堆積，抑制磷酸二酯酶的活性來(lái)刺激脂肪的分解，能有助于保護(hù)細(xì)胞免受紫外線輻射[39]；谷氨酸和谷氨酰胺等氨基酸是揮發(fā)性風(fēng)味物質(zhì)的前體物質(zhì)，影響肉質(zhì)風(fēng)味[40-41]。因此，這些差異代謝物可能是導(dǎo)致蜀興1號(hào)肉兔和伊拉兔肉質(zhì)差異的重要原因。

本研究通過(guò)對(duì)蜀興1號(hào)肉兔和伊拉兔的肉質(zhì)特性的比較分析，發(fā)現(xiàn)蜀興1號(hào)肉兔肌纖維密度和脂肪含量更高，并通過(guò)轉(zhuǎn)錄組學(xué)和代謝組學(xué)分析，篩選到SMTNL1、PM20D2和EDN1等可能影響肉質(zhì)的差異表達(dá)基因以及2-十一烯醛、肉堿、谷氨酸等差異代謝物，為蜀興1號(hào)肉兔的優(yōu)質(zhì)養(yǎng)殖以及優(yōu)質(zhì)肉兔新品種培育提供參考和依據(jù)。

參考文獻(xiàn)：

[1] 武拉平，王建勛，秦應(yīng)和. 2022年兔產(chǎn)業(yè)生產(chǎn)概況、2023年發(fā)展趨勢(shì)及政策建議[J]. 中國(guó)畜牧雜志，2023，59（3）：348-352.

[2] 謝躍杰. 伊拉兔肉腥味物質(zhì)及其體內(nèi)沉積變化的研究[D]. 重慶：西南大學(xué)，2016.

[3] 李叢艷，梅秀麗，鄺良德，等. 限飼對(duì)蜀興1號(hào)商品兔生長(zhǎng)性能、屠宰性能及肉品質(zhì)的影響[J]. 中國(guó)畜牧雜志，2022，58（1）：247-251.

[4] MENG X， GAO Z， LIANG Y， et al. Longissimus dorsi muscle transcriptomic analysis of Simmental and Chinese native cattle differing in meat quality[J]. Frontiers in Veterinary Science，2020，7：601064.

[5] CHEN B， YUE Y， LI J， et al. Transcriptome-metabolome analysis reveals how sires affect meat quality in hybrid sheep populations[J]. Frontiers in Nutrition，2022，9：967985.

[6] HU Z， CAO J， GE L， et al. Characterization and comparative transcriptomic analysis of skeletal muscle in Pekin duck at different growth stages using RNA-Seq[J]. Animals，2021，11（3）：834.

[7] REN L， LIU A， WANG Q， et al. Transcriptome analysis of embryonic muscle development in Chengkou Mountain Chicken[J]. BMC Genomics，2021，22（1）：431.

[8] OTT K H， ARANIBAR N， SINGH B， et al. Metabonomics classifies pathways affected by bioactive compounds. Artificial neural network classification of NMR spectra of plant extracts[J]. Phytochemistry，2003，62（6）：971-985.

[9] ZHOU H， YANG Y， WANG L， et al. Integrated multi-omic data reveal the potential molecular mechanisms of the nutrition and flavor in Liancheng white duck meat[J]. Frontiers in Genetics，2022，13：939585.

[10]CARTONI MANCINELLI A， SILLETTI E， MATTIOLI S， et al. Fatty acid profile， oxidative status， and content of volatile organic compounds in raw and cooked meat of different chicken strains[J]. Poultry Science，2021，100（2）：1273-1282.

[11]CHEN C， ZHENG J， XIONG C， et al. Metabolomics characterize the differential metabolic markers between Bama Xiang pig and Debao pig to identify pork[J]. Foods，2022，12（1）：5.

[12]劉 浪，熊國(guó)遠(yuǎn)，朱秀柏. 家兔的胴體性狀和肉品質(zhì)測(cè)定方法[J]. 中國(guó)養(yǎng)兔雜志，2009（3）：11-14.

[13]PERTEA M， PERTEA G M， ANTONESCU C M， et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads[J]. Nature Biotechnology，2015，33（3）：290-295.

[14]LIAO Y， SMYTH G K， SHI W. FeatureCounts：an efficient general purpose program for assigning sequence reads to genomic features[J]. Bioinformatics，2014，30（7）：923-930.

[15]ROBINSON M D， MCCARTHY D J， SMYTH G K. EdgeR：a Bioconductor package for differential expression analysis of digital gene expression data[J]. Bioinformatics，2010，26（1）：139-140.

[16]KIM D， LANGMEAD B， SALZBERG S L. HISAT：a fast spliced aligner with low memory requirements[J]. Nature Methods，2015，12（4）：357-360.

[17]李金玲，白昊鑫，馬晶晶，等. 肉鵝制品中松香酸和脫氫松香酸殘留的高效液相色譜檢測(cè)方法[J]. 江蘇農(nóng)業(yè)科學(xué)，2023，51（11）：174-178.

[18]梁 鵬，張 穩(wěn)，馮登偵，等. 基于轉(zhuǎn)錄組學(xué)篩選綿羊肉質(zhì)性狀相關(guān)候選基因[J]. 華北農(nóng)學(xué)報(bào)，2022，37（4）：220-231.

[19]陳 勇，朱友軍，鄒萬(wàn)瑞，等. 舍飼牦牛、犏牛與放牧牦牛的肌纖維特性研究[J]. 湖北農(nóng)業(yè)科學(xué)，2019，58（21）：141-144，149.

[20]PAN P， QIN Z， XIE W， et al. Identification of differentially expressed genes in the longissimus dorsi muscle of Luchuan and Duroc Pigs by transcriptome sequencing[J]. Genes （Basel），2023，14（1）：132.

[21]ZHAO L， LI F， ZHANG X， et al. Integrative analysis of transcriptomics and proteomics of longissimus thoracis of the Hu sheep compared with the Dorper sheep[J]. Meat Science，2022，193：108930.

[22]趙旺生，李柯銳，張婷婷，等. 基于高通量轉(zhuǎn)錄組測(cè)序的牦牛和犏牛附睪尾部差異表達(dá)基因分析[J]. 南方農(nóng)業(yè)學(xué)報(bào)，2023，54（5）：1273-1282.

[23]丁玥竹，任 宇，扈孟雪，等. 云南獨(dú)龍雞和紅原雞卵巢組織轉(zhuǎn)錄組水平的比較分析[J]. 江蘇農(nóng)業(yè)科學(xué)，2023，51（5）：38-44.

[24]袁 建，敖 政，曾素梅，等. 基于轉(zhuǎn)錄組測(cè)序的摘除卵巢藏姜母豬背最長(zhǎng)肌脂代謝相關(guān)基因篩選與分析[J]. 南方農(nóng)業(yè)學(xué)報(bào)，2023，54（3）：669-680.

[25]王錢保，姜潤(rùn)深，黎壽豐，等. 光照調(diào)控優(yōu)質(zhì)肉雞性成熟的轉(zhuǎn)錄組分析[J]. 江蘇農(nóng)業(yè)科學(xué)，2022，50（3）：16-21.

[26]熊 訊，阮 涌，許厚強(qiáng). 基于轉(zhuǎn)錄組分析干擾FABP1基因?qū)ωi肌內(nèi)脂肪沉積的影響[J]. 南方農(nóng)業(yè)學(xué)報(bào)，2023，54（3）：724-734.

[27]BORMAN M A， FREED T A， HAYSTEAD T A， et al. The role of the calponin homology domain of smoothelin-like "（SMTNL1） in myosin phosphatase inhibition and smooth muscle contraction[J]. Molecular and Cellular Biochemistry，2009，327（1/2）：93-100.

[28]LONTAY B， BODOOR K， SIPOS A， et al. Pregnancy and smoothelin-like protein "（smtnl1） deletion promote the switching of skeletal muscle to a glycolytic phenotype in human and mice[J]. Journal of Biological Chemistry，2015，290（29）：17985-17998.

[29]MURALI M， MACDONALD J A. Smoothelins and the control of muscle contractility[J]. Advances in Pharmacological Sciences，2018，81：39-78.

[30]MAJOR E， GYǒRY F， HORVTH D， et al. Smoothelin-like protein "regulates development and metabolic transformation of skeletal muscle in hyperthyroidism[J]. Frontiers in Endocrinology，2021，12：751488.

[31]HUANG X， HE P， WU L. Clinical significance of peptidase m20 domain containing "ii patients with carotid atherosclerosis[J]. Arquivos Brasileiros de Cardiologia，2022，119（3）：372-379.

[32]VEIGA-DA-CUNHA M， CHEVALIER N， STROOBANT V， et al. Metabolite proofreading in carnosine and homocarnosine synthesis：molecular identification of PM20D2 as β-alanyl-lysine dipeptidase[J]. Journal of Biological Chemistry，2014，289（28）：19726-19736.

[33]KOBAYASHI Y， YOSHIMOTO Y， YAMAMOTO Y， et al. Roles of EDNs in regulating oviductal NO synthesis and smooth muscle motility in cows[J]. Reproduction，2016，151（6）：615-622.

[34]張 潤(rùn)，楊 曼，王立賢，等. 畜禽肉中代謝物質(zhì)對(duì)肉品質(zhì)的影響及相關(guān)基因研究進(jìn)展[J]. 畜牧獸醫(yī)學(xué)報(bào)，2022，53（8）：2444-2452.

[35]SHI X， ZHANG X， SONG S， et al. Identification of characteristic flavour precursors from enzymatic hydrolysis-mild thermal oxidation tallow by descriptive sensory analysis and gas chromatography-olfactometry and partial least squares regression[J]. Journal of Chromatography B，2013，913/914：69-76.

[36]WANG Y R， WANG S L， LUO R M. Evaluation of key aroma compounds and protein secondary structure in the roasted Tan mutton during the traditional charcoal process[J]. Frontiers in Nutrition，2022，9：1003126.

[37]KAFFARNIK S， PREUSS S， VETTER W. Direct determination of flavor relevant and further branched-chain fatty acids from sheep subcutaneous adipose tissue by gas chromatography with mass spectrometry[J]. Journal of Chromatography A，2014，1350：92-101.

[38]DURAZZO A， LUCARINI M， NAZHAND A， et al. The nutraceutical value of carnitine and its use in dietary supplements[J]. Molecules，2020，25（9）：2127.

[39]HERMAN A， HERMAN A P. Caffeine’s mechanisms of action and its cosmetic use[J]. Skin Pharmacology and Physiology，2013，26（1）：8-14.

[40]LI Y， WANG J， WANG T， et al. Differences between Kazak Cheeses fermented by single and mixed strains using untargeted metabolomics[J]. Foods，2022，11（7）：966.

[41]WENG K， SONG L， BAO Q， et al. Comparative characterization of key volatile compounds in slow- and fast-growing duck raw meat based on widely targeted metabolomics[J]. Foods，2022，11（24）：3975.

（責(zé)任編輯：石春林）