基于RNA-Seq的向日葵脂肪酸合成基因篩選與分析

郭樹春,張艷芳,李素萍,邵 盈,邢慧軍,于海峰*

(1 內蒙古自治區農牧業科學院 作物科學研究所,呼和浩特 010031;2 內蒙古農業大學 園藝與植物保護學院,呼和浩特 010010;3 四子王旗林業和草原局,內蒙古烏蘭察布 011800)

植物油脂是由脂肪酸和甘油化合而成的天然高分子化合物。它是人體必需營養成分,可為人體提供必要的能量、維持特定的功能,同時它還是脂溶性維生素的載體[1]。脂肪酸中的單不飽和油酸、多不飽和亞油酸都是有益脂肪酸,其含量及組成是決定其營養價值和用途的關鍵[2]。向日葵(HelianthusannuusL.)是世界四大油料作物之一。向日葵籽油(又稱葵花籽油)是日常消費、食品加工中的一種優質植物油,其脂肪酸(fatty acid,FA)主要由油酸、亞油酸、棕櫚酸和硬脂酸組成[3]。向日葵籽油飽和脂肪酸(棕櫚酸和硬脂酸)相對含量很低,介于4%～13%之間,不飽和脂肪酸(油酸、亞油酸)含量高,在87%～94%之間[4-5],不飽和脂肪酸中油酸和亞油酸含量呈負相關關系。早在1999年聯合國糧農組織根據向日葵油酸含量的不同,將其基因型分為三大類:①低油酸向日葵(傳統向日葵),油酸含量14%～39%;②中油酸向日葵,油酸含量42%～72%;③高油酸向日葵,油酸含量75%～91%。現在市場上的向日葵油原料主要是傳統向日葵,其不飽和脂肪酸中亞油酸占55%～65%,油酸占20%～30%[6]。油酸因其分子結構中不飽和烯鍵數目少于亞油酸而更加穩定,高油酸向日葵油烹制的食品貨架期會更長,更加健康[7-8],因此市場對高油酸向日葵油更加青睞。現已證實,向日葵籽油中油酸和亞油酸的占比是決定其品質的關鍵,其脂肪酸含量和組分由基因型決定,還受到一些環境因子如溫度、光照和濕度等影響[9-10],且這些環境因子中對其影響最大的是溫度[11-12]。進一步研究表明,隨著溫度的升高,低油酸向日葵籽中的油酸含量明顯升高,亞油酸含量隨之降低,而溫度越低,則相反。而高油酸向日葵中的油酸含量對環境變化的敏感性較低,其脂肪酸組成變化也較小[10]。但是,從轉錄組層面在不同生長條件下對向日葵脂肪酸代謝的分子機制研究較少。

隨著中國農業的高質量發展,不同脂肪酸組分表型的向日葵需求將不斷增加。為了迎合社會發展的需要,彌補相關研究的不足,本研究以高油酸自交系J9、低油酸自交系P50為試驗材料,采用錯期播種方法,通過轉錄組分析,解析向日葵脂肪酸合成及組分變化的分子調控機制,促進高油酸向日葵分子育種穩步發展,為向日葵產業高質量發展提供理論依據。

1 材料和方法

1.1 試驗材料

以自育高油酸向日葵自交系J9(油酸相對含量大于85.0%)和低油酸自交系P50(油酸含量小于35.0%)為試驗材料,于2022年4月20日開始分期播種,每間隔10 d設1個播期,共設6個不同的播種期,每個播期處理設3次重復。試驗在呼和浩特市(N40.49°,E111.41°,海拔106 m)內蒙古自治區農牧業科學院試驗田進行,試材采取相同的栽培管理措施。

1.2 樣品采集

在開花授粉后(days after flowering,DAF)20 d,取整齊一致單株的向日葵籽仁為試驗樣品,累計采集12份樣品,液氮速凍。編號為:J9_S1、J9_S2、J9_S3、J9_S4、J9_S5、J9_S6和P50_S1、P50_S2、P50_S3、P50_S4、P50_S5、P50_S6,每份樣品采集6份,3份(3次重復)用于轉錄組測序,3份(3次重復)用于脂肪酸含量與組分測定。采樣期氣象數據來源于中國氣象數據網(http://data.cma.cn/)。

1.3 脂肪酸組分與含量測定

利用液相色譜法測定各樣品(累計36個)主要脂肪酸相對含量,如:油酸、亞油酸、硬脂酸、亞麻酸、棕櫚酸等組分參照GB 5009.168—2016第三法測定其含量。

1.4 轉錄組測序

利用TaKaRa公司的RNAisoPlus提取樣本總RNA[13],并文庫構建,構建好的文庫用Illumina HiSeqTM 4000進行轉錄組測序。

1.5 測序數據過濾與組裝

使用Trinity軟件對轉錄組的De novo進行從頭組裝[14]。原始序列數據經雜質去除后得到clean reads,然后將每個樣品clean reads用短reads比對軟件TMAP比對到參考基因序列上(DDBJ/ENA/GenBank,MNCJ00000000[15]),組裝得到的Unigene[16],得到最終的Unigene信息,具體信息分析過程及分析流程見圖1。

圖1 數據處理流程

1.6 基因表達量計算和功能注釋

使用RPKM法(reads per kb per million reads)[17]計算基因的表達量,兩樣本間的差異表達基因篩選參考Audic[18]基于測序的差異基因檢測方法,此方法發表于GenomeResearch,公式為:MRPKM=(1000000×C)/(N×L/1000)。通過控制FDR(fa1se discovery rate)≤0.001且︱log2R︱≥1(倍數差異在2倍以上)的值校正P閾值。

RPKM法公式各符號含義:A基因的表達量用RPKM(A)表示,則C表示能唯一比對到基因A區域的Reads數,N為唯一比對到參考基因的總reads數,L為基因A的堿基數。具體計算過程:首先將測序得到的reads匹配到參考基因上,將能匹配到參考基因上的reads數標準化到106,以此為基礎來換算某一基因上匹配上的reads數,再除以該基因長度就得到該基因的表達量(RPKM值)。RPKM不僅對測序深度作了歸一化,而且對基因長度也作了歸一化,使得不同長度的基因在不同測序深度下得到的基因表達水平估計值有了可比性,是目前最可靠的基因表達估計方法,計算得到的基因表達量可直接用于比較不同樣品間的基因表達差異分析。

每個樣品中的每個Unigene表達量確定后,為了便于后續對基因的功能研究,對每個基因進行功能注釋,將每個基因序列首先跟NCBI的NR庫用(http://www.geneontology.org/)BLAST軟件進行比對,再將每個基因注釋到NT、Swiss-Prot、KEGG、COG、GO功能庫。例如:用Blast2 GO軟件注釋到GO的各級條目下。基因序列通過BLAST比對到KEGG數據庫進行KEGG注釋,Q≤0.05的Pathway定義為在差異表達基因中顯著富集的Pathway,通過Pathway顯著性富集能確定差異表達基因參與的主要生物代謝途徑和信號轉導通路[19]。

2 結果與分析

2.1 油酸含量測定結果

錯期播種同期采樣試驗條件下,測定自交系J9和P50在6個播期的油酸相對含量,并收集采樣期(從開始授粉到20 DAF)平均溫度。結果(圖2)顯示:高油酸品系J9油酸相對含量受日平均溫度變化影響不顯著,其油酸相對含量基本穩定在85.0%左右;而低油酸品系P50油酸相對含量隨日平均溫度的降低而顯著降低,在S1播種期油酸相對含量最高達42.0%,而S6播種期油酸相對含量僅為20.49%。由此可見,向日葵脂肪酸組分和含量易受溫度影響,且高油酸向日葵品系脂肪酸組分和含量受溫度影響小,而低油酸品系易受溫度影響。因此,筆者設定了錯期播種同期采樣試驗,在授粉期不同溫度下篩選向日葵脂肪酸合成基因,并研究其調控關系。

*表示兩品系在同一時期的顯著水平(P<0.05)。下同。

2.2 RNA-Seq測序質量評估結果

通對36個(12個樣品,3個重復)樣品使用Hiseq 4000測序平臺進行RNA-Seq測序,2個品系共測得265.79 G數據,平均每個樣品數據約7.38 G,單個樣品數據最少的為J9_S2_2,共測得6.12 G數據;最多的為P50_S2_2,共測得9.06 G數據。樣品平均Q20為96.59%,平均Q30為92.93%。堿基錯誤率分布符合一般規律。GC含量平均值為46.50%,且不存在明顯分離現象,結果見表1。上述結果顯示,RNA-Seq測序數據質量達到預期,可用于后續分析。

表1 測序數據質量評估

2.3 RNA-Seq測序重復性檢驗

根據每個樣本基因表達量的數值,計算皮爾遜相關系數的平方(R2),得出重復樣本之間的相關性。結果(圖3)為:在J9品系中,相關系數最高的為J9_S5_1和J9_S5_3之間,相關系數高達0.935;相關系數最低的為J9_S1_2和J9_S1_3之間,相關系數為0.754。在P50品系中,相關系數最高的為P50_S4_1和P50_S4_3之間,相關系數高達0.918;相關系數最低的為P50_S3_2和P50_S3_1之間,相關系數為0.769。平均相關系數高于0.8。因此,本研究中各個生物學重復間相關性較好,樣品的重復性良好。

圖中橫向、縱向數字1,2,3表示3個重復。

2.4 Unigene基因表達量分析和功能注釋

利用測序得到的clean reads,采用RPKM算法確定基因的表達量。首先將某一樣品測序得到的reads匹配到參考基因上,并將reads數標準化為106,以此為基礎來換算該基因上匹配上的reads數,再除以該基因長度得到該樣品該基因的基因表達量,每個樣本每個基因的基因表達量用RPKM方法表示。按RPKM值的大小,將基因分為高表達基因(RPKM≥10)、中表達基因(1

圖4 不同向日葵品系基因表達量的疊加直方圖

對上述所有Unigene進行功能注釋,這些基因被注釋到了GO功能庫下的1 283個GO條目下。在3個GO本體(ontology)中,生物學過程(biological process,BP)中條目最多,包含665個條目;其次分子功能(molecular function,MF)和細胞組分(cellular component,CC),分別包含451,166個條目。在圖5中列出了位于3個GO本體中前10位的GO條目,KEGG庫注釋的Unigene共計3 286個,其中位于前20位中的角質層、木栓質和蠟生物合成(cutin, suberine and wax biosynthesis)屬于脂質代謝。

圖5 所有Unigene GO和Unigene KEGG分類結果

2.5 差異表達基因篩選與分析結果

首先對品種內不同處理樣品差異表達分析,將同一品種不同播期品與其對照(S1)間的比對所獲得的差異表達基因總和定義為S-CK數據集,S-CK數據集(S2-CK、S3-CK、S4-CK、S5-CK、S6-CK,累計5個數據集,表2)反映了J9和P50在脂肪酸積累關鍵時期,在不同播期下條件下,品種內響應環境變化的DEGs;把J9和P50在相同播期的測序結果經比對所獲得的差異表達基因之和定義為S數據集(累計6個數據集,表2),該數據集反映J9和P50在脂肪酸積累各關鍵時期,品種間響應環境變化的DEGs。

表2 差異表達基因篩選數據集命名

對S-CK數據集進行品種內差異分析發現,在S-CK數據集中,高油酸品系J9在S5-CK中的差異基因數(different expression genes,DEGs)最多,為1 022個;低油酸品系P50在S3-vs-S1中的差異基因數最多,為1 712個。上述結果表明,高油酸品系在S5響應環境變換的差異表達基因較多,而低油酸品系在S3期響應環境變換的差異表達基因較多;總體上高油酸品系響應環境變換的差異表達基因較少,這和2.1節向日葵籽粒中油酸含量和溫度變化關系研究結果相吻合。

在基因表達調控類型研究中,S-CK數據集的5種比較(S2 vs S1、S3 vs S1、S4 vs S1、S5 vs S1和S6 vs S1)中,高油酸品系J9中含有258,382,100,311,207個上調DEGs,有551,467,346,711,1 427個下調DEGs;低油酸品系P50中含有12,930,34,305,873個上調DEGs,有12,782,42,355,980個下調DEGs(圖6,A和圖6,B);上述與油酸合成相關DEGs的獲得,為油酸合成的共表達分析以及通路分析提供了基礎。

A.品種內上調差異分析結果;B.品種內下調差異分析結果;C.品種內特異表達分析。粉色填充表上調;綠色填充表下調;黃色線表J9品系;藍色線表P50品系。

上述結果經試驗內特異表達分析顯示,在S-CK數據集的5種比較(S2 vs S1、S3 vs S1、S4 vs S1、S5 vs S1和S6 vs S1)策略中,高油酸品系J9中含有372,279,115,322,147個特異上調DEGs;有664,386,718,660,1 249個特異下調DEGs(圖6,C),分析發現J9和P50的品種特異性表達基因數遠超其品種共表達數量,該結論證明了本研究材料選擇的可靠性,提取特異表達基因功能條目,為后續研究提供信息。同時,為了進一步確定影響高油酸和低油酸品系的關鍵差異基因,對S數據集(來源于RNA-Seq數據)進行品種間特異表達分析(見表2),結果如圖7。由圖7可知,在S數據集中,無論是上調基因還是下調基因,各時期品種間特異表達基因相當。最后,品種特異性表達基因和品種間特異表達基因去除重復后,在J9和P50在-CK和S數據集得到15 885個與向日葵脂肪酸代謝相關的品種間特異表達DEGs用于后續功能分析。

圖7 J9 vs P50在相同播種期差異表達基因分析

2.6 特異表達基因功能分析

將2.5節分析得到的15 885個特異表達DEGs映射到GO注釋中,經GO功能富集,得到與脂質代謝相關的GO條目:如,脂質生物合成過程(GO:0008610)、脂肪酸代謝過程(GO:0006631)、類固醇生物合成過程(GO:0006694)和類固醇代謝過程(GO:0008202)等。J9在中特異表達DEGs的GO注釋中,累計獲得96與脂質代謝相關DEGs(表3),P50中獲得43個(表4);J9和P50在品種間特異表達GO條目映射研究中,獲得278個DEGs(表5)。

表3 J9特異表達DEGs GO term映射結果

表4 P50特異表達DEGs GO term映射結果

表5 J9和P50特異表達DEGs GO term映射結果

總之,在品種內特異表達分析和品種間特異表達分析中,累計403個DEGs被注釋到與脂質代謝相關的GO條目下。在品種內特異表達分析結果中,將通過GO功能富集獲得的與脂質代謝直接相關的403個DEGs映射到與脂質代謝相關的KEGG代謝通路,累計29個DEGs直接參與了19條相關KEGG通路(圖8,A)。其中,顯著富集于糖基磷脂酰肌醇(GPI)生物合成KEGG通路的DEGs 5個、脂肪酸延伸的4個、脂肪酸代謝的4個、類固醇生物合成的3個、脂肪酸生物合成的3個、萜類骨架生物合成的3個、不飽和脂肪酸生物合成的2個、α-亞麻酸代謝的2個和甘油磷脂代謝通路的3個(圖8)。

圖8 特異表達DEGs的KEGG功能富集

3 討 論

將在品種間特異表達分析結果與品種內特異表達分析得到DEGs和去除重復DEGs,共得到42個DEGs,經分析其隸屬于22種酶。以此為基礎繪制了向日葵脂肪酸代謝相關DEGs通路圖,向日葵脂肪酸代謝過程論述如下。

首先,向日葵通過光合作用合成蔗糖(suc),然后被蔗糖轉用酶(suc transporter,SUC)運送到向日葵種子細胞中,經過磷酸化、糖酵解、氧化生成脂肪酸合成的碳源乙酰輔酶A(acetyl-CoA),在質體(plastid)中開始其脂肪酸的從頭合成(圖9)。首先,乙酰-CoA在乙酰-CoA羧化酶(ACCase)催化下,與酰基載體蛋白(acyl carrier protein,ACP)結合,生成向日葵脂肪酸合成的直接原料丙二酰-CoA(malonyl-CoA);然后,丙二酰-CoA在脂肪酸合成酶復合體(fatty acid synthase,FAS)或酰基載體蛋白轉移酶(MCAT)催化下[20]生成脂肪酸碳鏈延長供體丙二酰-ACP(malonyl-ACP);同時,乙酰-CoA,ACP轉酰基酶將1個丙二酰基團由乙酰-CoA轉移至ACP的絲氨酸殘基上,生成乙酰-ACP,并釋放CoA(圖10)。本試驗中,基因110889119和110925594參與了脂肪酸合成的限速酶乙酰-CoA羧化酶(ACCase)的調節,基因110889119在J9和P50的脂肪酸積累前期相對表達量偏高,其在J9中更高;基因110925594P50卻在P50中的相對表達量更高,說明該基因更容易受到環境因素的影響。

圖9 脂肪酸生物合成通路中的表達模式

其次,由β-酮酯酰-ACP還原酶(β-ketoacyl-ACP reductase,KAR)催化乙酰-ACP生成β-ketoacyl-ACP,然后在NADPH等酶的參與下,生成β-羥基丁酰-ACP(β-hydroxyacyl-ACP),此步是脂肪酸生物合成的最初還原步驟。本試驗中,有累計6個基因參與了上述過程的調節,其中調節β-酮酯酰-ACP還原酶(KAR)的2個基因在A-試驗中表達特征明顯,基因110864403在積累前期低表達,后期高表達,而基因110893771則表現相反。4個基因參與了NADPH等酶的調節,基因1109382561和110868833表達特征明顯,基因110868833在整個播期在J9中表達量較高,而1109382561則相反,它們在A-試驗中也表現出同樣的特征;值得注意的是,雖然基因110868833在J9整個脂肪酸積累期表達量相對較高,但是其在油酸積累關鍵期(20 DAFs)表達量最高。

再次,β-羥基丁酰-ACP由β-羥基丁酰-ACP脫水酶(β-hydroxacyl-ACP dehydratase,HAD)催化[21],在其α位和β位碳原子之間脫水生成2,3-反式丁烯酰(巴豆酰)-ACP(enoyl-ACP)。然后,反式丁烯酰-ACP被多功能蛋白(multifunctional protein,MFP)催化[22],轉化生成相應的丁酰-ACP(butyryl-ACP)。基因110893751參與了此過程中MFP的激活,該基因在J9中脂肪酸合成28 DAFs和35 DAFs表達量最高。

最后,飽和丁酰-ACP再經轉酰基作用與新的丙二酸單酰-ACP依次重復上述4個過程,經數次循環后,脂肪酸合成通過終止反應來終止合成,可生成C16∶0和C18∶1飽和脂肪酸;然后,在脂酰基-ACP硫脂酶(acyl-ACP thioesterase,Fat)作用下,將ACP上的酰基部分(C16∶0和C18∶1飽和脂肪酸)通過水解反應釋放出來,成為游離脂肪酸,并釋放1個巰基-ACP,使反應終止[23]。碳鏈聚合終止后,產生的不同碳鏈長度游離脂肪酸,由酰基-CoA合成酶(acyl-CoA sythetase,ACS)負責催化合成脂酰-CoA(C16∶0-CoA和C18∶1-CoA),并采用某種未知機制從質體轉運到內質網或胞質中,完成脂肪酸鏈的延長。此過程中,基因110890752和110871748參與了C16∶0和C18∶1飽和脂肪酸的合成,而基因110941185和110893751參與了脂酰-CoA的合成與轉運。

4 結 論

在向日葵籽粒形成過程中,雖然其脂肪酸組分和含量隨著積累時間的推移而變化,但是向日葵脂肪酸主要以不飽和脂肪酸油酸和亞油酸為主,其合成的關鍵轉折時期是20 DAF;通過播種期調控可以明顯影響向日葵脂肪酸組分和含量,向日葵脂肪酸組分和含量易受溫度影響,且高油酸向日葵材料脂肪酸組分和含量受溫度影響小,而低油酸材料易受溫度影響。在錯期播種試驗條件下,對36份樣品進行轉錄組測序,2個品系共測得265.79 G數據。

試驗中,J9在S5-vs-S1中的差異基因數最多,為1 022個;P50在S3-vs-S1中的差異基因數最多,為1 712個。經基因差異表達分析,得到的15 885個特異表達DEGs,經GO功能富集,得到與脂質代謝相關的GO條目累計獲得96與脂質代謝相關DEGs;最后,J9和P50品種間特異表達GO條目映射研究中,獲得278個DEGs。

為進一步揭示向日葵種子發育過程中與油脂積累及油酸生物合成相關的關鍵基因,將通過GO功能富集獲得的與脂質代謝直接相關的403個DEGs映射到與脂質代謝相關的KEGG代謝通路,累計29個DEGs直接參與了19條相關KEGG通路;去除重復DEGs,共得到43個DEGs,隸屬于22種酶。其中KAS、ACX、DGAT、MGD和OLE是向日葵不飽和脂肪酸生物合成和油脂積累的關鍵基因。