不同品種紫花苜蓿轉錄組分析及營養品質差異的探討

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

(1.內蒙古農業大學草原與資源環境學院，農業部飼草栽培、加工與高效利用重點試驗室，草地資源教育部重點實驗室，內蒙古 呼和浩特 010011；2.內蒙古自治區科技信息研究院，內蒙古 呼和浩特 010010；3.內蒙古自治區草原勘察規劃院，內蒙古 呼和浩特 010051；4.內蒙古自治區農牧業科學院，內蒙古 呼和浩特 010031；5.四川省會東縣農牧局，四川 涼山 615000)

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

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

在本研究中，使用Illumina測序技術，在初花期對種植在同一塊試驗地的2個紫花苜蓿品種(國外品種：WL319HQ；國內品種：準格爾)的葉片RNA樣品進行denovo組裝，得到紫花苜蓿的Unigenes數據，從而獲得整個紫花苜蓿葉片組織的轉錄組信息。通過分析其差異基因和代謝通路，試圖從基因層面解釋造成不同品種營養差異的內部機理。通過該項研究為紫花苜蓿生產實踐和遺傳育種奠定理論基礎，同時為后續的紫花苜蓿轉錄組研究提供參考。

1 材料與方法

1.1 試驗地概況和試驗材料

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

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

1.2 實驗流程和測定方法

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

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

1.2.3denovo組裝和Unigenes注釋 由于以前沒有獲得紫花苜蓿基因組信息，因此使用參考基因組獨立的Trinity方法[11]進行denovo組裝得到Unigenes。使用BLASTx程序(http://www.ncbi.nlm.nih.gov/BLAST/)將Unigenes序列比對到蛋白數據庫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。根據最佳比對結果獲得蛋白質功能注釋。

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

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

1.3 數據處理方法

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

2 結果與分析

2.1 不同品種營養品質比較

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

2.2 測序和組裝

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

表1 2個苜蓿品種葉片營養品質比較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 不同苜蓿品種樣品測序數據統計結果Table 2 Statistics of sequencing data of the different alfalfa varieties samples

2.3 de novo組裝結果的評估分析

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

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

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

2.4 基因注釋結果分析

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

2.5 差異基因的篩選

表5 基因注釋結果統計Table 5 Statistics of gene annotation results

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

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

圖2 2個紫花苜蓿葉片的差異基因分析火山圖 Fig.2 Volcano map of differential gene statistics of 2 alfalfa leaves 橫坐標表示兩個樣品的差異倍數對數值，縱坐標表示兩個樣品的FDR的負log10值。紅色(WL319HQ相對于準格爾表達量上調)和綠色(表達量下調)的點表示基因的表達量有差異(判斷標準為 FDR<0.05, 差異倍數兩倍以上)，黑色的點為沒有差異。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)中，“代謝過程”，“細胞過程”，“刺激反應”和“單一生物過程”富集基因相對較多，分別有288，260，174和159個Unigenes。其中“代謝過程”中有195個上調Unigenes和93個下調Unigenes，“細胞過程”中有170個上調Unigenes和90個下調Unigenes，“刺激反應”中有131個上調Unigenes和43個下調Unigenes，“單一生物過程”中的Unigenes有109個上調，50個下調。而“免疫系統過程”和“節律過程”富集的基因較少，只有2個Unigenes。在細胞組分(cellular component)中，大多數基因被注釋到“細胞”(210個上調，46個下調)，“細胞部分”(210個上調，46個下調)和“細胞器”(156個上調，30個下調)。對于分子功能(molecular function)，“結合”(169個上調，97個下調)和“催化活性”(124個上調，99個下調)是富集基因較多的功能類別。

2.7 KEGG途徑富集分析

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

表6 GO功能統計表Table 6 GO functional statistics table

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

3 討論

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

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

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

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

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

4 結論

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