成啟明，格根圖，撒多文，王志軍，范文強(qiáng)，卜振鯤，司強(qiáng)，李俊峰，盧娟，賈玉山*

(1.內(nèi)蒙古農(nóng)業(yè)大學(xué)草原與資源環(huán)境學(xué)院，農(nóng)業(yè)部飼草栽培、加工與高效利用重點(diǎn)試驗(yàn)室，草地資源教育部重點(diǎn)實(shí)驗(yàn)室，內(nèi)蒙古 呼和浩特 010011；2.內(nèi)蒙古自治區(qū)科技信息研究院，內(nèi)蒙古 呼和浩特 010010；3.內(nèi)蒙古自治區(qū)草原勘察規(guī)劃院，內(nèi)蒙古 呼和浩特 010051；4.內(nèi)蒙古自治區(qū)農(nóng)牧業(yè)科學(xué)院，內(nèi)蒙古 呼和浩特 010031；5.四川省會(huì)東縣農(nóng)牧局，四川 涼山 615000)

紫花苜蓿(Medicagosativa)作為世界上應(yīng)用最廣、經(jīng)濟(jì)價(jià)值最高和栽培面積最大的多年生優(yōu)質(zhì)豆科牧草，在現(xiàn)代草業(yè)的可持續(xù)發(fā)展中具有突出的地位[1]。紫花苜蓿莖葉中含有豐富的蛋白質(zhì)、礦物質(zhì)、多種氨基酸、維生素和胡蘿卜素等，其總能、消化能、代謝能也較高，被認(rèn)為是高蛋白優(yōu)質(zhì)飼草[2]。作為全球性的飼料作物，苜蓿產(chǎn)業(yè)的發(fā)展在一定程度上反映著一個(gè)國家或地區(qū)畜牧業(yè)發(fā)展的水平和質(zhì)量。因此，對苜蓿進(jìn)行全面的營養(yǎng)價(jià)值評價(jià)是畜牧業(yè)發(fā)展的必要條件。紫花苜蓿品種多種多樣，受自身遺傳特性及環(huán)境條件的限制，其適宜生長的環(huán)境和長勢也就有所不同，其內(nèi)部的營養(yǎng)特性自然也會(huì)存在一定的差異。趙燕梅等[3]通過對11個(gè)苜蓿品種的營養(yǎng)特性分析發(fā)現(xiàn)，不同品種間的干物質(zhì)、粗蛋白、中性洗滌纖維、酸性洗滌纖維、可溶性糖等營養(yǎng)物質(zhì)差異極顯著。于輝等[4]在哈爾濱地區(qū)通過對4個(gè)紫花苜蓿品種進(jìn)行大田試驗(yàn)發(fā)現(xiàn)，阿爾岡金苜蓿的草產(chǎn)量高于其他品種，和平苜蓿的營養(yǎng)品質(zhì)較好，肇東苜蓿的越冬率最高，適合在該地區(qū)種植?？悼∶返萚5]在北京地區(qū)對10個(gè)紫花苜蓿品種進(jìn)行長達(dá)6年的引種試驗(yàn)發(fā)現(xiàn)，愛菲尼特、牧歌和勝利者3個(gè)苜蓿品種的年均干草產(chǎn)量顯著高于其他品種，適合在北京地區(qū)種植。不同品種苜蓿不但在產(chǎn)量和營養(yǎng)物質(zhì)含量上存在很大差異，高微微等[6]通過對45個(gè)紫花苜蓿品種的代謝產(chǎn)物皂苷進(jìn)行測定研究發(fā)現(xiàn)，不同品種間皂苷含量差異顯著。而目前關(guān)于從其遺傳特性出發(fā)，從基因?qū)用嫣接懺斐善錉I養(yǎng)品質(zhì)差異的內(nèi)部機(jī)理的研究還鮮有報(bào)道。

所有表達(dá)基因的身份以及其轉(zhuǎn)錄水平，綜合起來被稱作轉(zhuǎn)錄組(Transcriptome)，廣義上指某一生理?xiàng)l件下細(xì)胞內(nèi)所有轉(zhuǎn)錄產(chǎn)物的集合，包括信使RNA、核糖體RNA、轉(zhuǎn)運(yùn)RNA及非編碼RNA；狹義上指所有信使RNA的集合[7]。轉(zhuǎn)錄組測序研究在新基因的發(fā)現(xiàn)，基因功能注釋，基因差異表達(dá)和分子標(biāo)記發(fā)展中具有重要價(jià)值[8-10]。通過新一代Illumina測序技術(shù)，能夠全面快速地獲得某一物種特定組織或器官在某一狀態(tài)下的幾乎所有轉(zhuǎn)錄本序列信息，并已經(jīng)廣泛用于模式和非模式植物的研究。

在本研究中，使用Illumina測序技術(shù)，在初花期對種植在同一塊試驗(yàn)地的2個(gè)紫花苜蓿品種(國外品種：WL319HQ；國內(nèi)品種：準(zhǔn)格爾)的葉片RNA樣品進(jìn)行denovo組裝，得到紫花苜蓿的Unigenes數(shù)據(jù)，從而獲得整個(gè)紫花苜蓿葉片組織的轉(zhuǎn)錄組信息。通過分析其差異基因和代謝通路，試圖從基因?qū)用娼忉屧斐刹煌贩N營養(yǎng)差異的內(nèi)部機(jī)理。通過該項(xiàng)研究為紫花苜蓿生產(chǎn)實(shí)踐和遺傳育種奠定理論基礎(chǔ)，同時(shí)為后續(xù)的紫花苜蓿轉(zhuǎn)錄組研究提供參考。

1 材料與方法

1.1 試驗(yàn)地概況和試驗(yàn)材料

試驗(yàn)地點(diǎn)位于中國內(nèi)蒙古包頭市。地跨東經(jīng)10°37″-110°27″，北緯40°5″-40°17″。屬北溫帶大陸氣候，干旱多風(fēng)，春季干旱少雨多風(fēng)，夏季溫和短促，秋季涼爽溫差大，冬季溫長而寒冷。年平均氣溫6.8 ℃，7月平均氣溫22.5～23.1 ℃，1月平均氣溫-13.7 ℃。無霜期約165 d，最大凍土深度1.4 m。年平均降水量330 mm，年平均蒸發(fā)量2094 mm，日平均風(fēng)速3 m·s-1；全年日照時(shí)數(shù)3177 h，年日照百分率70%，是全國日照最豐富的地區(qū)之一。

試驗(yàn)材料為種植在同一試驗(yàn)地的兩個(gè)品系四倍體紫花苜蓿，中國品種(準(zhǔn)格爾)和美國品種(WL319HQ)，于2017年6月1日，取2個(gè)品種的初花期葉片，并按照說明書使用TRIzol(Invitrogen)分離總的RNA樣品。分別設(shè)置3個(gè)生物學(xué)重復(fù)，共計(jì)6個(gè)轉(zhuǎn)錄組樣品，本試驗(yàn)材料無參考基因組。

1.2 實(shí)驗(yàn)流程和測定方法

1.2.1RNA提取和RNA-Seq文庫構(gòu)建 使用RNeasy Plant Mini Kit(Qiagen，Germany)試劑盒按照制造商的說明書從2個(gè)紫花苜蓿品種的葉片中提取總RNA。通過Agilent 2100生物分析儀(Agilent Technologies，USA)檢測RNA樣品的質(zhì)量。通過NEBNext Oligo(dT)25珠(NEB，USA)從50 μL總RNA中富集Poly(A)mRNA。然后按照說明書，通過用于Illumina的NEBNext Ultra RNA文庫制備試劑盒(NEB，USA)將富集的mRNA構(gòu)建成cDNA文庫。

1.2.2Illumina測序 將提取的總RNA樣品直接送往廣州基迪奧生物科技有限公司使用Illumina HiSeq 4000測序平臺進(jìn)行轉(zhuǎn)錄組測序。測序得到的原始圖像數(shù)據(jù)先將base calling轉(zhuǎn)化為序列數(shù)據(jù)，同時(shí)過濾數(shù)據(jù)中這些雜質(zhì)Raw reads，從而得到Clean reads。

1.2.3denovo組裝和Unigenes注釋 由于以前沒有獲得紫花苜?；蚪M信息，因此使用參考基因組獨(dú)立的Trinity方法[11]進(jìn)行denovo組裝得到Unigenes。使用BLASTx程序(http://www.ncbi.nlm.nih.gov/BLAST/)將Unigenes序列比對到蛋白數(shù)據(jù)庫Nr (http://www.ncbi.nlm.nih.gov)、Swissprot (http://www.expasy.ch/sprot)、KEGG(http:// www.genome.jp/kegg)和KOG (http://www.genome.jp/kegg/ko.html)，其E值閾值為1×10-5。根據(jù)最佳比對結(jié)果獲得蛋白質(zhì)功能注釋。

1.2.4差異表達(dá)基因的選擇 測序后的基因使用Reads per kilobase per million mapped reads(RPKM)法計(jì)算基因表達(dá)量[12]。使用edgeR軟件包(http://www.r-project.org/)來選擇品種間的差異基因。差異表達(dá)倍數(shù)≥2倍的基因被認(rèn)定是差異基因[13]。對篩選出的差異表達(dá)基因進(jìn)行GO功能富集分析和KEGG途徑分析。根據(jù)Nr注釋信息，使用BLAST2GO[14]將苜蓿葉片轉(zhuǎn)錄本進(jìn)行功能分類，得到Unigene的GO注釋信息。得到每個(gè)Unigene的GO注釋后，用WEGO軟件[15]對所有Unigene做GO功能分類統(tǒng)計(jì)?；贙EGG途徑分析有助于更進(jìn)一步了解基因的生物學(xué)功能。

1.2.5營養(yǎng)指標(biāo)測定方法 粗蛋白(crude protein，CP)、中性洗滌纖維(neutral detergent fibre，NDF)、酸性洗滌纖維(acid detergent fiber，ADF)按照《飼料分析及飼料質(zhì)量檢測技術(shù)》測定[16]。

1.3 數(shù)據(jù)處理方法

利用 Microsoft Office Excel 2007 軟件進(jìn)行圖、表和數(shù)據(jù)的前期處理，利用SAS 9.1.3(statistical analysis system)軟件進(jìn)行數(shù)據(jù)計(jì)算及方差分析。

2 結(jié)果與分析

2.1 不同品種營養(yǎng)品質(zhì)比較

如表1所示，苜蓿品種對各營養(yǎng)指標(biāo)含量的影響較大，WL319HQ苜蓿的粗蛋白(CP)含量顯著高于準(zhǔn)格爾苜蓿(P<0.05)，而中性洗滌纖維(NDF)和酸性洗滌纖維(ADF)含量顯著低于準(zhǔn)格爾苜蓿(P<0.05)。

2.2 測序和組裝

如表2所示，為了比較2個(gè)紫花苜蓿品種的差異基因，取2個(gè)紫花苜蓿品種的初花期葉片，產(chǎn)生了6個(gè)DGE文庫，具有3個(gè)生物學(xué)重復(fù)。使用HiSeq 4000對所有文庫進(jìn)行測序，得到約39G的總的核苷酸， 2億多個(gè)reads。 其中每個(gè)DGE文庫得到約6.5G的核苷酸，約4300萬個(gè)reads，其中質(zhì)量不低于20的堿基的比例(Q20)均超過了96%，不能測序的核苷酸“N”的百分比為0，堿基G和C數(shù)占總堿基數(shù)的百分比大于42%。通過Trinity組裝得到66734個(gè)Unigenes，平均長度為869 bp。

表1 2個(gè)苜蓿品種葉片營養(yǎng)品質(zhì)比較Table 1 Comparison of nutritional quality of leaves of two alfalfa varieties (DM%)

注：同列不同小寫字母表示差異顯著(P<0.05)。

Note： Different lowercase letters within the same column indicate significant different at 0.05 level.

表2 不同苜蓿品種樣品測序數(shù)據(jù)統(tǒng)計(jì)結(jié)果Table 2 Statistics of sequencing data of the different alfalfa varieties samples

2.3 de novo組裝結(jié)果的評估分析

如表3所示，總共有26474，24230和28904個(gè)Unigenes分別匹配到鷹嘴豆(Cicerarietinum)，大豆(Glycinemax)，蒺藜苜蓿(Medicagotruncatula)和紫花苜蓿的直系同源編碼序列。其中有590(2.23%)，457(1.89%)和343(1.19%)個(gè)Unigenes的長度分別與鷹嘴豆，大豆和蒺藜苜蓿的同源基因相等(ratio=1)；有3001(11.34%)，3317(13.69%)和2785(9.64%)個(gè)Unigenes的長度分別小于鷹嘴豆，大豆和蒺藜苜蓿的同源基因(ratio<1)；有22883(86.43%)，20456(84.42%)和25776(89.17%)個(gè)Unigenes的長度分別大于鷹嘴豆，大豆和蒺藜苜蓿的同源基因(ratio>1)。如表4所示，9636(58.71%)個(gè)鷹嘴豆同源基因，6382(34.68%)個(gè)大豆同源基因和8195(49.06%)個(gè)蒺藜苜蓿同源基因可以被Unigenes覆蓋，覆蓋率大于80%；覆蓋百分比為50%～80%的3367(20.52%)個(gè)鷹嘴豆同源基因，4810(26.13%)個(gè)大豆同源基因和2969(17.77%)個(gè)蒺藜苜蓿同源基因被Unigenes覆蓋；另外有2085(12.71%)個(gè)鷹嘴豆同源基因，3779(20.53%)個(gè)大豆同源基因和3330(19.93%)個(gè)蒺藜苜蓿同源基因能夠被Unigenes覆蓋，覆蓋百分比為20%～50%；此外有1322(8.06%)個(gè)鷹嘴豆同源基因，3434(18.66%)個(gè)大豆同源基因和2785(13.24%)個(gè)蒺藜苜蓿同源基因僅被20%的Unigenes覆蓋。

表3 Unigene長度/直系同源物長度的比率Table 3 Ratio of Unigene length/ortholog length

表4 Unigene/ortholog覆蓋百分比Table 4 Cover percentage of Unigene/ortholog

2.4 基因注釋結(jié)果分析

如表5所示，總共有45657(68.42%)條Unigenes被注釋到，而有21077(31.58%)條Unigenes沒有被注釋。對于Nr數(shù)據(jù)庫，通過BLASTx注釋到44888(67.26%)條Unigenes，對于Swissprot，KOG和KEGG數(shù)據(jù)庫分別注釋到29190(43.74%)，24844(37.23%)和15647(23.45%)條Unigenes。有11690條Unigenes被同時(shí)注釋到4個(gè)蛋白數(shù)據(jù)庫，有47條Unigenes只注釋到KEGG數(shù)據(jù)庫，40條Unigenes單獨(dú)注釋到KOG數(shù)據(jù)庫，12156條Unigenes單獨(dú)注釋到Nr數(shù)據(jù)庫，有287條Unigenes單獨(dú)注釋到Swissprot數(shù)據(jù)庫(圖1)。

2.5 差異基因的篩選

表5 基因注釋結(jié)果統(tǒng)計(jì)Table 5 Statistics of gene annotation results

通過對2個(gè)不同紫花苜蓿品種的基因測序， 在錯(cuò)誤發(fā)現(xiàn)率(false discovery rate，F(xiàn)DR)<0.05 且 |log2 fold change (FC)|>1的篩選條件下找到2個(gè)紫花苜蓿品種的差異表達(dá)基因。如圖2所示，在2個(gè)紫花苜蓿葉片中找到1098個(gè)差異基因 (DEGs) ，其中有706個(gè)Unigenes上調(diào)，392個(gè)Unigenes下調(diào)(準(zhǔn)格爾-vs-WL319HQ)。得到的差異表達(dá)基因用于后續(xù)分析。

圖1 由BLASTx注釋的Unigenes數(shù)量的維恩圖Fig.1 Venn diagram of number of Unigenes annotated by BLASTx

圖2 2個(gè)紫花苜蓿葉片的差異基因分析火山圖 Fig.2 Volcano map of differential gene statistics of 2 alfalfa leaves 橫坐標(biāo)表示兩個(gè)樣品的差異倍數(shù)對數(shù)值，縱坐標(biāo)表示兩個(gè)樣品的FDR的負(fù)log10值。紅色(WL319HQ相對于準(zhǔn)格爾表達(dá)量上調(diào))和綠色(表達(dá)量下調(diào))的點(diǎn)表示基因的表達(dá)量有差異(判斷標(biāo)準(zhǔn)為 FDR<0.05, 差異倍數(shù)兩倍以上)，黑色的點(diǎn)為沒有差異。The abscissa represents the logarithm of the difference multiples of the two samples, and the ordinate represents the negative log10 value of FDR of the two samples. The points of red (WL319HQ is up-regulated relative to Zhungeer) and green (down-regulated) indicate that the expression levels of the genes are difference (the standard is FDR<0.05, the difference multiple is more than twice), while the black points showed no difference.

2.6 GO功能顯著性富集分析

如表6所示，通過GO功能注釋，在生物過程(biological process)中，“代謝過程”，“細(xì)胞過程”，“刺激反應(yīng)”和“單一生物過程”富集基因相對較多，分別有288，260，174和159個(gè)Unigenes。其中“代謝過程”中有195個(gè)上調(diào)Unigenes和93個(gè)下調(diào)Unigenes，“細(xì)胞過程”中有170個(gè)上調(diào)Unigenes和90個(gè)下調(diào)Unigenes，“刺激反應(yīng)”中有131個(gè)上調(diào)Unigenes和43個(gè)下調(diào)Unigenes，“單一生物過程”中的Unigenes有109個(gè)上調(diào)，50個(gè)下調(diào)。而“免疫系統(tǒng)過程”和“節(jié)律過程”富集的基因較少，只有2個(gè)Unigenes。在細(xì)胞組分(cellular component)中，大多數(shù)基因被注釋到“細(xì)胞”(210個(gè)上調(diào)，46個(gè)下調(diào))，“細(xì)胞部分”(210個(gè)上調(diào)，46個(gè)下調(diào))和“細(xì)胞器”(156個(gè)上調(diào)，30個(gè)下調(diào))。對于分子功能(molecular function)，“結(jié)合”(169個(gè)上調(diào)，97個(gè)下調(diào))和“催化活性”(124個(gè)上調(diào)，99個(gè)下調(diào))是富集基因較多的功能類別。

2.7 KEGG途徑富集分析

如表7所示，有8414個(gè)基因注釋到87個(gè)代謝通路中， 這些基因中有281個(gè)差異基因。其中Unigenes數(shù)量注釋排名前10的代謝通路分別是：“核糖體”(ko03010)，“剪接體”(ko03040)，“內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工”(ko04141)，“碳代謝”(ko01200)，“內(nèi)吞作用”(ko04141)，“氨基酸的生物合成”(ko01230)，“RNA轉(zhuǎn)運(yùn)”(ko03013)，“植物-病原互作”(ko04626)，“淀粉和蔗糖代謝”(ko00500)，“植物激素信號轉(zhuǎn)導(dǎo)”(ko04075)。其中“核糖體”途徑注釋的DEGs最多，為64 (22.78%)，依次為“內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工”(44 個(gè)DEGs，占15.66%)、“碳代謝”(22個(gè)DEGs，占7.83%)和“剪接體”(19 個(gè)DEGs，占6.76%)。

表6 GO功能統(tǒng)計(jì)表Table 6 GO functional statistics table

表7 Unigenes 數(shù)量排名前10的代謝通路Table 7 Top 10 metabolic pathways involving Unigenes

3 討論

苜蓿作為國內(nèi)外優(yōu)質(zhì)蛋白質(zhì)飼草，其適口性和營養(yǎng)價(jià)值都很高。而苜蓿原料中粗蛋白(CP)、中性洗滌纖維(NDF)和酸性洗滌纖維(ADF)的含量多少是反映其品質(zhì)好壞的重要指標(biāo)，CP含量越高，NDF和ADF含量越低，其飼草品質(zhì)就越好[17]。在本研究中WL319HQ苜蓿的CP含量顯著高于準(zhǔn)格爾苜蓿，而NDF和ADF含量顯著低于準(zhǔn)格爾苜蓿，說明WL319HQ的品質(zhì)優(yōu)于準(zhǔn)格爾。本研究發(fā)現(xiàn)，不同品種紫花苜蓿的營養(yǎng)品質(zhì)差異顯著，該研究結(jié)果與趙燕梅等[3]的研究結(jié)果一致，造成其營養(yǎng)品質(zhì)差異可能是因?yàn)椴煌贩N的遺傳特性差異所導(dǎo)致的，而目前關(guān)于從基因和代謝通路層面出發(fā)探討不同苜蓿品種營養(yǎng)品質(zhì)差異的研究還鮮有報(bào)道。因此本研究開展了不同品種紫花苜蓿的轉(zhuǎn)錄組測序，想要從基因和代謝通路層面探討其營養(yǎng)品質(zhì)差異的內(nèi)在機(jī)理。

轉(zhuǎn)錄組測序是在高通量測序技術(shù)基礎(chǔ)上發(fā)展起來的轉(zhuǎn)錄組分析技術(shù)，其中應(yīng)用比較廣泛的轉(zhuǎn)錄組測序平臺是Illumines測序平臺。Illumina測序技術(shù)為研究基因結(jié)構(gòu)和功能提供了強(qiáng)大的技術(shù)支持，特別是對于無參考基因組的非模式植物進(jìn)行的denovo從頭測序，可以獲得該物種的參考序列，從而為其后續(xù)的基因研究提供理論基礎(chǔ)。目前關(guān)于利用Illumina法對沒有參考基因組物種進(jìn)行轉(zhuǎn)錄組測序的研究較多，其中包括速生桉(Eucalyptusrobustasmith)[18]、甘薯(Dioscoreaesculenta)[19]、胡蘿卜(Daucuscarota)[20]和胡黃連(Picrorhizascrophulariiflora)[21]等。本試驗(yàn)利用高通量測序平臺Illumina HiSeq 4000對2個(gè)紫花苜蓿品種(無參考基因組)葉片進(jìn)行轉(zhuǎn)錄組分析。基于denovo組裝，總共獲得66734條Unigenes，平均長度為869 bp，本研究得到的Unigenes數(shù)量明顯高于劉希強(qiáng)等[22]關(guān)于紫花苜蓿轉(zhuǎn)錄組測序的研究結(jié)果(41734條Unigenes)，而低于張森浩[23]的研究結(jié)果(192875條Unigenes)。造成這種差異的原因可能是品種和生長環(huán)境的不同。本研究獲得的紫花苜蓿Unigenes平均長度明顯高于鵝掌楸(Liriodendronchinensis)(537 bp)[24]、短絲木犀(Osmanthusserrulatus)(697 bp)[25]和桂花(Osmanthus)(708 bp)[26]等，與火力楠(Micheliamacclurei)(852 bp)[27]和青葉膽(Herbaswertiaemileensis)(901 bp)[28]等相近。總體來講，本研究得到的轉(zhuǎn)錄本和Unigenes長度與已有的紫花苜蓿轉(zhuǎn)錄組測序研究結(jié)果相似[29-30]，說明本研究測序和拼接結(jié)果符合預(yù)期，其轉(zhuǎn)錄組數(shù)據(jù)可信。在本研究中有些Unigenes不能完全覆蓋鷹嘴豆，大豆和蒺藜苜蓿的同源基因的完整編碼序列，但是本研究中2個(gè)紫花苜蓿品種的Unigenes能夠覆蓋大部分的同源基因，并且大部分Unigenes長度與鷹嘴豆，大豆和蒺藜苜蓿同源基因長度的比值約等于1，說明本研究的轉(zhuǎn)錄組序列質(zhì)量較好。

將66734條Unigenes序列通過使用BLASTx(評估<10-5)，對Nr，Swissprot，KOG和KEGG的4個(gè)公開蛋白質(zhì)數(shù)據(jù)庫進(jìn)行搜索。從注釋的結(jié)果中可以發(fā)現(xiàn)，45657(68.42%)個(gè)Unigenes至少在一個(gè)數(shù)據(jù)庫得到注釋，大約有11690個(gè)Unigenes被同時(shí)注釋到4個(gè)蛋白數(shù)據(jù)庫。其中，Nr數(shù)據(jù)庫注釋到的Unigenes數(shù)量最多，為44888(67.26%)條，這與鄧敏捷等[31]對泡桐樹(Paulowniasieb)進(jìn)行轉(zhuǎn)錄組測序的研究結(jié)果(65.5%)相似。然而，在本研究中產(chǎn)生了21077條沒有被注釋的Unigenes，很可能代表著參試種的一個(gè)特殊基因群。

目前基于轉(zhuǎn)錄組測序平臺來研究不同品種紫花苜蓿差異基因的報(bào)道很少。本研究利用轉(zhuǎn)錄組測序找到了2個(gè)紫花苜蓿品種的1098個(gè)差異基因 (DEGs)，相對于國內(nèi)品種準(zhǔn)格爾而言，國外苜蓿WL319HQ苜蓿中有706個(gè)上調(diào)基因，392個(gè)下調(diào)基因?；谛蛄邢嗨菩裕瑢?098個(gè)DEGs分配給生物過程(biological process)、細(xì)胞組分(cellular component)和分子功能(molecular function)的一個(gè)或多個(gè)本體。在生物過程中，“代謝過程”(288個(gè)Unigenes)和“細(xì)胞過程”(260個(gè)Unigenes)富集的Unigenes占主體地位，這與張雪松等[32]的研究結(jié)果一致。而“免疫系統(tǒng)過程”(2個(gè)Unigenes)和“節(jié)律過程”(2個(gè)Unigenes)富集的基因較少。在細(xì)胞組分中，大多數(shù)基因被注釋到“細(xì)胞”(256個(gè)Unigenes)，“細(xì)胞部分”(256個(gè)Unigenes)和“細(xì)胞器”(186個(gè)Unigenes)。對于分子功能，“結(jié)合”(266個(gè)Unigenes)和“催化活性”(223個(gè)Unigenes)是富集基因較多的功能類別。

KEGG數(shù)據(jù)庫是全面研究基因產(chǎn)物在細(xì)胞中的代謝途徑和研究基因產(chǎn)物的功能的數(shù)據(jù)庫，基于KEGG的代謝途徑的分析有助于更進(jìn)一步了解基因的生物學(xué)功能[32]。在本研究中，有8414個(gè)基因注釋到87個(gè)代謝通路中，這些基因中有281個(gè)差異基因。其中Unigenes 數(shù)量注釋較多的途徑是“核糖體”(622個(gè)Unigenes，64個(gè)DEGs)，“剪接體”(562個(gè)Unigenes，19個(gè)DEGs)，“內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工”(550個(gè)Unigenes，44個(gè)DEGs)，“碳代謝”(484個(gè)Unigenes，22個(gè)DEGs)。粗蛋白質(zhì)(CP)是常規(guī)飼料分析中用以估計(jì)飼料、動(dòng)物組織或動(dòng)物排泄物中一切含氮物質(zhì)的指標(biāo)，它包括了真蛋白質(zhì)和非蛋白質(zhì)含氮物(non-protein nitrogen, NPN)兩部分。NPN包括游離氨基酸、硝酸鹽、氨等。蛋白質(zhì)由20種氨基酸組成，其生物合成可分為5個(gè)階段：氨基酸的活化、多肽鏈合成的起始、肽鏈的延長、肽鏈的終止和釋放、蛋白質(zhì)合成后的加工修飾。其合成場所為核糖體，在本研究中“核糖體”代謝途徑的注釋基因最多，其差異基因也是最多，這可能是導(dǎo)致參試品種蛋白差異的主要原因之一。而“剪接體”、“內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工”和“碳代謝”途徑也是直接或者間接影響蛋白質(zhì)合成的代謝途徑，因此，有可能是這些代謝途徑的差異，導(dǎo)致這2個(gè)品種紫花苜蓿的蛋白質(zhì)含量差異。不同的組織細(xì)胞具有不同的生理功能，是因?yàn)樗鼈儽磉_(dá)不同的基因，產(chǎn)生具有特殊功能的蛋白質(zhì)，參與蛋白質(zhì)生物合成的成份至少有200種，其主要體是由mRNA、tRNA、核糖核蛋白體以及有關(guān)的酶和蛋白質(zhì)因子共同組成。因此蛋白質(zhì)的合成與轉(zhuǎn)化是一個(gè)極其復(fù)雜的過程，是多個(gè)代謝途徑和多個(gè)基因聯(lián)合作用的，因此關(guān)于不同品種苜蓿粗蛋白質(zhì)含量差異的原因有待更加深入的研究。在本研究中，有294個(gè)基因注釋到“植物激素信號轉(zhuǎn)導(dǎo)”(ko04075)代謝途徑，其中有6個(gè)差異基因。而“植物激素信號轉(zhuǎn)導(dǎo)”代謝途徑主要是控制植物激素的合成和轉(zhuǎn)化。植物激素在植物的生長發(fā)育過程中起到至關(guān)重要的作用，植物的生長發(fā)育伴隨著細(xì)胞的分裂和分化，從而導(dǎo)致植物細(xì)胞壁組分和含量也伴隨著變化[33]。植物細(xì)胞壁的主要組成成分是纖維素、半纖維素、木質(zhì)素和果膠。因此，“植物激素信號轉(zhuǎn)導(dǎo)”代謝途徑間接地影響著中性洗滌纖維(NDF)和酸性洗滌纖維(ADF)的合成，這可能是導(dǎo)致本研究中2個(gè)紫花苜蓿品種的NDF和ADF含量差異的重要原因之一。本研究從基因和代謝通路層面初步探討出導(dǎo)致2個(gè)苜蓿品種營養(yǎng)品質(zhì)差異的內(nèi)在原因，然而具體的營養(yǎng)指標(biāo)差異的內(nèi)在機(jī)理有待進(jìn)一步深入研究和驗(yàn)證。

4 結(jié)論

本研究通過Illumina HiSeq 4000平臺，首次利用RNA-seq技術(shù)研究了不同紫花苜蓿品種(準(zhǔn)格爾、WL319HQ)的差異表達(dá)基因，初步建立了它們的基因表達(dá)譜。在2個(gè)紫花苜蓿品種中找到1098個(gè)差異基因 (DEGs) ，其中有706個(gè)Unigenes上調(diào)，392個(gè)Unigenes下調(diào)(準(zhǔn)格爾-vs-WL319HQ)。對其差異基因進(jìn)行GO功能分析和KEGG途徑分析，初步分析了2個(gè)品種營養(yǎng)品質(zhì)差異的內(nèi)在原因。本研究極大地豐富了紫花苜蓿轉(zhuǎn)錄組數(shù)據(jù)信息，為今后紫花苜蓿的轉(zhuǎn)錄組測序提供理論依據(jù)，同時(shí)為紫花苜蓿生產(chǎn)實(shí)踐提供參考。