整合關(guān)聯(lián)分析和共表達(dá)網(wǎng)絡(luò)分析挖掘甘藍(lán)型油菜籽粒質(zhì)量候選基因

陳元軍, 馬娟娟, 史 睿, 李偉龍, 王 婷, 彭 琦, 張 維, 陳 鋒, 王曉東, 高建芹, 付三雄, 張潔夫, 孫程明, 季彪俊, 胡茂龍,3,4

(1.福建農(nóng)林大學(xué)農(nóng)學(xué)院,福建 福州 350002; 2.江蘇省農(nóng)業(yè)科學(xué)院經(jīng)濟(jì)作物研究所/農(nóng)業(yè)農(nóng)村部長江下游棉花與油菜重點實驗室/江蘇省農(nóng)業(yè)生物學(xué)重點實驗室/江蘇省現(xiàn)代作物生產(chǎn)協(xié)同創(chuàng)新中心,江蘇 南京 210014; 3.南京農(nóng)業(yè)大學(xué)作物遺傳與種質(zhì)創(chuàng)新國家重點實驗室,江蘇 南京 210095; 4.江蘇大學(xué)生命科學(xué)研究院,江蘇 鎮(zhèn)江 212023)

甘藍(lán)型油菜(BrassicanapusL.)是中國主要的食用油來源之一,種植面積約7.0×106hm2,所產(chǎn)植物油量約占國產(chǎn)植物油總量的55%。但中國菜籽油產(chǎn)量已連續(xù)十年低于消費(fèi)量,2020年,菜籽油產(chǎn)量為6.36×106t,消費(fèi)量為8.25×106t,缺口達(dá)1.89×106t,這對中國糧食安全構(gòu)成了重大威脅[1]。因此,油菜產(chǎn)量的提高成為當(dāng)今油菜育種工作中最為關(guān)鍵的目標(biāo)。籽粒質(zhì)量是油菜產(chǎn)量性狀的重要組成部分,相較于角果長、每角粒數(shù)等其他產(chǎn)量性狀具有更高的遺傳力[2]。研究結(jié)果表明,籽粒質(zhì)量改良是近年來油菜產(chǎn)量提升的重要原因[3]。因此,解析籽粒質(zhì)量遺傳基礎(chǔ)、挖掘有利基因可以提高油菜產(chǎn)量、保障中國油料供給安全。

甘藍(lán)型油菜是雙子葉植物,油菜種子由胚、胚乳和種皮等部分組成,分別由受精卵、受精中央細(xì)胞和珠被組織發(fā)育而成。基因型通過調(diào)控胚乳發(fā)育決定種子大小,珠被組織決定了種子發(fā)育的空間,受精卵發(fā)育成胚,三者結(jié)合,共同決定了種子的質(zhì)量和大小[4]。目前種子發(fā)育大小的主要形成機(jī)制分為2類,一類是種子早期發(fā)育階段胚乳增殖形成的空腔大小決定了種子體積,隨后胚乳細(xì)胞分裂促進(jìn)種子質(zhì)量增加;另一類是種皮和角果皮等母體組織調(diào)控種子發(fā)育,種子在發(fā)育過程中依賴母體的光合作用提供的養(yǎng)分[5]。目前,一系列基因已被報道調(diào)控種子發(fā)育,如IKU1通過調(diào)控胚乳生長影響種子籽粒質(zhì)量[6];STK通過控制種皮的結(jié)構(gòu)和機(jī)械性能調(diào)節(jié)細(xì)胞分裂素水平來控制果實大小[7];GW5通過泛素蛋白酶體途徑調(diào)節(jié)種子細(xì)胞分裂,控制種子粒寬和籽粒質(zhì)量等[8]。

甘藍(lán)型油菜籽粒質(zhì)量是由多基因控制的數(shù)量性狀,受基因型和環(huán)境共同作用[9]。王瑞等[10]利用7個黃籽油菜品系配制完全雙列雜交組合,利用ADM(加性-顯性-母體遺傳效應(yīng))遺傳模型分析產(chǎn)量性狀,結(jié)果表明,千粒質(zhì)量的遺傳主要由基因的加性和顯性效應(yīng)聯(lián)合決定,而母體效應(yīng)比重小,同時環(huán)境對各性狀的遺傳均有一定影響。崔嘉成等[11]利用11個甘藍(lán)型油菜親本雜交得到30個組合,對各組合及親本角果性狀的遺傳模型的分析結(jié)果表明,基因的加性效應(yīng)控制千粒質(zhì)量的遺傳,無顯性效應(yīng)。戚存扣等[12]對2個遺傳差異大的親本構(gòu)建的世代家系群體的分析結(jié)果表明,該組合千粒質(zhì)量主要受主基因的加性效應(yīng)和多基因的加性和顯性效應(yīng)共同決定。李娜[13]的研究結(jié)果表明,油菜籽粒質(zhì)量主要受母體基因型調(diào)控,籽粒質(zhì)量的自然變異由細(xì)胞數(shù)目的變異決定,而非由種子密度決定。盡管油菜籽粒質(zhì)量的遺傳基礎(chǔ)十分復(fù)雜,但前人研究形成的共識是籽粒質(zhì)量遺傳力較高,籽粒質(zhì)量主要受到基因間加性效應(yīng)的影響,顯性效應(yīng)影響較小,雜種優(yōu)勢不強(qiáng)。因此,可以通過聚合不同種質(zhì)中的有利基因來增加油菜籽粒質(zhì)量。

全基因組關(guān)聯(lián)分析(Genome-wide association study,GWAS)是一種以位點間連鎖不平衡為基礎(chǔ),用統(tǒng)計學(xué)分析檢測遺傳多態(tài)性與目標(biāo)性狀間關(guān)聯(lián)的分析方法[14]。如今GWAS已被廣泛應(yīng)用于解析油菜的重要產(chǎn)量性狀。Lu等[15]對520份甘藍(lán)型油菜種質(zhì)資源的7個產(chǎn)量性狀進(jìn)行關(guān)聯(lián)分析,共檢測到128個顯著位點,其中93個為新檢測到的位點,結(jié)合4個產(chǎn)量差異極端材料的轉(zhuǎn)錄組數(shù)據(jù),確定14個候選基因。Dong等[16]對157份甘藍(lán)型油菜的籽粒質(zhì)量進(jìn)行關(guān)聯(lián)分析,檢測到20個籽粒質(zhì)量顯著位點,在A09染色體顯著位點附近找到油菜籽粒質(zhì)量和角果長已克隆基因BnaA.ARF18.a。Khan等[17]對521份甘藍(lán)型油菜品種的胚珠數(shù)、每角粒數(shù)和籽粒質(zhì)量性狀進(jìn)行關(guān)聯(lián)分析,利用多位點和單位點模型分別檢測到280個和31個顯著單核苷酸多態(tài)性(SNP),并挖掘出42個產(chǎn)量性狀候選基因。林升麗[18]通過重測序?qū)?05份油菜種質(zhì)資源進(jìn)行基因分型,對關(guān)聯(lián)群體7個環(huán)境的籽粒質(zhì)量表型進(jìn)行關(guān)聯(lián)分析,檢測到127個顯著位點,找到BnaA07.KCR、BnaA07.EXPA1和BnaA10.GIF2等多個候選基因。

加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(Weighted gene co-expression network analysis,WGCNA)是利用大量表達(dá)譜數(shù)據(jù)對表達(dá)模式類似的基因進(jìn)行聚類,挖掘與性狀關(guān)聯(lián)基因模塊的分析方法。WGCNA目前已被廣泛應(yīng)用于油菜重要性狀的候選基因挖掘。Ma等[19]基于油菜感染黑脛病病菌的轉(zhuǎn)錄組數(shù)據(jù),對檢測到的1 892個差異剪接基因進(jìn)行WGCNA,檢測到8個共表達(dá)基因模塊,并挖掘出屬于WRKY、AP2、HMG和C2H2等轉(zhuǎn)錄因子家族的樞紐基因。同時,越來越多的研究將GWAS與WGCNA相結(jié)合,找到一系列調(diào)控目標(biāo)性狀的候選基因。鮮小華等[20]利用GWAS挖掘出1 826個黃籽名義候選基因,利用WGCNA找到2個黃籽相關(guān)基因模塊,挖掘到BnF3H、BnANS等多個候選基因。劉景森[21]利用GWAS檢測到14個生物產(chǎn)量顯著位點,利用葉片、種子和莖稈等組織的轉(zhuǎn)錄組數(shù)據(jù)構(gòu)建油菜生物產(chǎn)量共表達(dá)網(wǎng)絡(luò),鑒定到BnA04g06420D、BnA09g35380、BnC03g73810D和BnC08g48810D等多個生物產(chǎn)量候選基因。田貴福等[22]整合了油菜籽粒質(zhì)量的GWAS和QTL定位結(jié)果,利用71份油菜的角果轉(zhuǎn)錄組數(shù)據(jù)構(gòu)建了共表達(dá)網(wǎng)絡(luò),挖掘出樞紐基因BnSWH1并驗證了其對種子籽粒質(zhì)量的負(fù)調(diào)控作用。本研究采用整合GWAS和WGCNA的策略,挖掘控制籽粒質(zhì)量的候選基因,為今后籽粒質(zhì)量相關(guān)的基因克隆和調(diào)控機(jī)制研究提供參考。

1 材料與方法

1.1 材料與數(shù)據(jù)來源

本研究的關(guān)聯(lián)群體包含496份國內(nèi)外的地方品種、育成品種及高世代育種材料,其中有444份材料來源于國內(nèi)油菜主產(chǎn)區(qū),包括江蘇、湖北、湖南和四川等,國外材料來自德國、瑞典和加拿大等國家。試驗材料于2014年種植于江蘇南京(記作14NJ),2015年、2016年種植于江蘇泰州(記作15TZ、16TZ)。在前期研究中,我們已利用MLM(混合線性模型)和GLM(一般線性模型)對該群體進(jìn)行籽粒質(zhì)量關(guān)聯(lián)分析,本研究所使用的基因型和表型源于前期發(fā)表的論文[23]。

本研究中用于WGCNA的轉(zhuǎn)錄組數(shù)據(jù)集來自國家基因組科學(xué)數(shù)據(jù)中心(National Genomics Data Center,NGDC)的生物項目數(shù)據(jù)庫BioProject(ID:PRJCA001246),包含大粒品種中雙11和中小粒品種中油821開花后第7 d、第10 d、第14 d、第45 d的種子和開花后第7 d、第10 d、第14 d的角果皮共計28份轉(zhuǎn)錄組數(shù)據(jù)。

1.2 全基因組關(guān)聯(lián)分析

本研究利用A-D Test模型對籽粒質(zhì)量基因型和表型進(jìn)行關(guān)聯(lián)分析,顯著性閾值(P)計算方法為0.01/總標(biāo)記數(shù),取對數(shù)為-lgP=6.28。若1 Mb區(qū)間內(nèi)存在較多的顯著SNP,把兩兩間決定系數(shù)(R2)≥0.1的SNP視為1個關(guān)聯(lián)位點,選取其中P值最小的為代表,利用qqman包繪制QQ圖(Quantile-quantile plot)和曼哈頓圖(Manhattan plot)[24]。

1.3 權(quán)重基因共表達(dá)網(wǎng)絡(luò)分析

利用WGCNA包構(gòu)建基因共表達(dá)網(wǎng)絡(luò)[25],過濾變異系數(shù)(CV)<0.4、表達(dá)量(TPM)≤1和缺失樣本比例超過10%的基因,根據(jù)基因間的表達(dá)量數(shù)據(jù)對樣品進(jìn)行層次聚類分析,采取逐步法聚類表達(dá)模式相近的基因,繪制出聚類樹和每個樣品的表達(dá)圖。為了使結(jié)果符合無尺度網(wǎng)絡(luò)分布[26],選取擬合曲線R2接近0.9時的權(quán)重參數(shù)β,power值為10。利用動態(tài)剪枝法識別共表達(dá)模塊,計算每個模塊的特征向量,合并距離較近的模塊并選擇相關(guān)系數(shù)絕對值大于0.75的模塊作為研究的目標(biāo)模塊。使用omicshare在線分析平臺(https://www.omicshare.com/tools/)對目標(biāo)模塊進(jìn)行GO富集,對富集結(jié)果進(jìn)行柱狀圖、氣泡圖和圈圖等圖表的繪制和展示。將|模塊身份(Module membership,KME)|≥0.90的基因定義為樞紐基因,利用Cytoscape軟件對籽粒質(zhì)量相關(guān)模塊進(jìn)行網(wǎng)絡(luò)可視化展示[27]。

1.4 候選基因挖掘

以R2=0.2時的連鎖不平衡(Linkage disequilibrium,LD)衰減距離作為GWAS顯著位點的置信區(qū)間,提取區(qū)間內(nèi)的基因編碼序列(CDS),與模式植物擬南芥的基因CDS進(jìn)行BLAST比對(E-value≤1×10-10),用相似比最大的擬南芥基因來標(biāo)注油菜基因。同理,注釋W(xué)GCNA中籽粒質(zhì)量相關(guān)模塊的基因。最終,將2個分析中找到的油菜籽粒質(zhì)量已知基因和擬南芥籽粒質(zhì)量已知基因的同源基因作為候選基因[28]。

2 結(jié)果與分析

2.1 全基因組關(guān)聯(lián)分析

利用A-D Test模型對籽粒質(zhì)量最佳線性無偏預(yù)測值(BLUP)進(jìn)行關(guān)聯(lián)分析,設(shè)置顯著性閾值為-lgP=6.28,共檢測到37個顯著SNP;合并存在連鎖不平衡的 SNP(R2≥0.1),得到19個顯著SNP,分布在A01、A02、A06～A09、C02、C03、C05、C06、C09等11條染色體上(圖1、表1)。利用一般線性模型計算, 19個位點聯(lián)合解釋34.1%的表型變異。顯著性最高位點Bn-scaff_17526_1-p1066214位于C09染色體1.49 Mb處,-lgP=13.06(表2);顯著性次高位點Bn-scaff_16804_1-p203519位于C02染色體6.98 Mb處,-lgP=11.24;其余位點的-lgP值范圍為6.47～9.67。與單個環(huán)境的試驗結(jié)果相比,14個位點在至少2個環(huán)境中被重復(fù)檢測到,其中,顯著性排名前4的位點Bn-scaff_17526_1-p1066214、Bn-scaff_16804_1-p203519、Bn-scaff_17119_1-p148890和Bn-A08-p20343735在3個環(huán)境(14NJ、15TZ及16TZ)中均被重復(fù)檢測到,證實了本研究結(jié)果的可重復(fù)性和可靠性(表2)。

表1 籽粒質(zhì)量顯著關(guān)聯(lián)位點

表2 各環(huán)境下A-D Test模型關(guān)聯(lián)位點

A:籽粒質(zhì)量A-D Test曼哈頓圖;B:籽粒質(zhì)量A-D Test QQ圖。a:A01;b:A02;c:A03;d:A04;e:A05;f:A06;g:A07;h:A08;i:A09;j:A10;k:C01;l:C02;m:C03;n:C04;o:C05;p:C06;q:C07。圖1 油菜籽粒質(zhì)量A-D Test全基因組關(guān)聯(lián)分析Fig.1 Genome-wide association study of seed weight in rapeseed by A-D Test model

比較前人GLM和MLM模型結(jié)果,8個位點與GLM重疊,2個位點與MLM重疊(表1)。位點Bn-scaff_17526_1-p1066214和Bn-A02-p6564861在3個模型中都被重復(fù)檢測到。不僅如此,A-D test還找到11個新位點,Bn-scaff_16804_1-p203519分布在C02染色體6.98 Mb處,-lgP=11.24,在3個環(huán)境中都被檢測到。Bn-A09-p29991443位于A09染色體27.82 Mb處,-lgP=8.39,在2個環(huán)境(15TZ、16TZ)中都被檢測到。位點Bn-A01-p3061223和位點Bn-A08-p4273091均在14NJ和16TZ中被檢測到,-lgP值分別為7.98和7.75,此外還有7個新位點位于各染色體上,分別在1～2個環(huán)境中被檢測到。合并3個模型的結(jié)果,得到71個顯著位點,11個分布在A07染色體上,9個分布在C06染色體上,8個分布在A03染色體上,其余染色體上分布1～7個位點。分析顯著位點的分布,57.7%的位點分布在A亞基因組,42.3%的位點分布在C亞基因組。利用一般線性模型計算,發(fā)現(xiàn)71個位點聯(lián)合解釋50.1%的表型變異。

2.2 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析

大粒品種中雙11的千粒質(zhì)量約為4.6 g,中小粒品種中油821的千粒質(zhì)量約為3.4 g。本研究下載了兩者開花后第7 d、10 d、14 d、45 d的種子和開花后第7 d、10 d、14 d的角果皮轉(zhuǎn)錄組數(shù)據(jù),利用R軟件中WGCNA包對28份轉(zhuǎn)錄組數(shù)據(jù)的100 919個基因進(jìn)行分析,過濾低表達(dá)(TPM≤1)和變異系數(shù)較低(CV<0.4)的基因后,保留了35 540個基因。通過層次聚類繪制樣本樹并構(gòu)建無尺度網(wǎng)絡(luò)(圖2A、圖2B),采取逐步法聚類表達(dá)模式相近的基因(圖2C),最終得到13個共表達(dá)模塊。不同模塊間包含的基因數(shù)目差別大,turquoise模塊含有數(shù)量最多的基因,包含7 176個基因;salmon模塊包含基因數(shù)目最少,僅包含108個基因;其余模塊內(nèi)基因數(shù)目范圍為185～6 189。根據(jù)模塊間相關(guān)性繪制相關(guān)性熱圖,purple與green,purple與pink,magenta與yellow,magenta與black等模塊間相關(guān)性較高(圖2D)。

2.3 篩選目標(biāo)模塊及GO富集分析

計算模塊特征向量與籽粒質(zhì)量表型的相關(guān)系數(shù),結(jié)果顯示,模塊系數(shù)絕對值范圍為0.004 9～0.820 0,以0.75為閾值,得到2個籽粒質(zhì)量關(guān)聯(lián)模塊,其中purple模塊系數(shù)最高(|r|>0.80,P<4×10-7),magenta模塊次之(|r|>0.76,P<2×10-6)(圖3A)。

A:模塊與籽粒質(zhì)量性狀關(guān)聯(lián)熱圖;B:purple模塊表達(dá)值熱圖與特征值條形圖;C:magenta模塊表達(dá)值熱圖與特征值條形圖。圖3 基因共表達(dá)模塊分析Fig.3 Gene co-expression module analysis

purple模塊包含1 109個基因,分布在19條染色體上,各染色體基因數(shù)目為25～92。比較2個品種,purple模塊基因在大粒品種中雙11種子和角果皮中表達(dá)量較低,在中小粒品種中油821中表達(dá)量較高(圖3B)。對模塊基因進(jìn)行GO富集分析,共富集到108個顯著GO 途徑(圖4A),包括苯丙氨酸氨基轉(zhuǎn)移酶活性、天冬氨酸轉(zhuǎn)氨酶活性、2-酮戊二酸轉(zhuǎn)氨酶活性、硫雙加氧酶活性和肉桂酸生物合成等(圖4B)。這些分子功能在控制種子的生長發(fā)育中有重要作用,如苯丙氨酸參與幼苗早期細(xì)胞分裂生長和種子發(fā)育等重要過程[29];天冬氨酸轉(zhuǎn)氨酶(AAT)是合成氨基酸的一種關(guān)鍵酶,Zhou等[30]研究發(fā)現(xiàn)天冬氨酸轉(zhuǎn)氨酶在水稻中過表達(dá)影響種子的氮代謝和氨基酸含量;硫雙加氧酶是擬南芥基因ETHE1編碼的一種對種子胚胎和胚乳發(fā)育至關(guān)重要的酶類[31];Baleroni等[32]研究發(fā)現(xiàn)肉桂酸衍生物對甘藍(lán)型油菜種子子葉的脂肪酸和總脂質(zhì)積累有積極作用。

A:紫色(purple)模塊GO富集圈圖;B:purple模塊GO富集分析中生物過程分類前20個GO富集途徑;C:洋紅色(magenta)模塊GO富集氣泡圖。D:magenta模塊前20個GO富集途徑。Z值表示某個基因集合在特定生物過程或功能方面的顯著性程度。圖4 目標(biāo)模塊GO富集分析Fig.4 Target module GO enrichment analysis

magenta模塊含有1 114個基因,分布在19條染色體上,各染色體基因數(shù)目為23～127。與purple模塊相反,magenta模塊基因在大粒品種中雙11中表達(dá)豐度顯著高于中油821(圖3C)。對模塊基因進(jìn)行GO富集分析,共富集到521個顯著GO term。圖4C和圖4D展示的是magenta模塊中生物學(xué)過程、分子功能和細(xì)胞組分三大類TOP 20 GO通路氣泡圖,其中核酸代謝過程、雜環(huán)代謝過程、細(xì)胞氮化合物代謝過程、細(xì)胞芳香族化合物代謝過程等代謝通路被顯著富集。這些通路的基因參與了種子的發(fā)育過程,如在模式植物擬南芥中,焦磷酸酶的表達(dá)促進(jìn)種子中蛋白質(zhì)和油類物質(zhì)的合成,UPF1RNA解旋酶參與調(diào)控種子大小[33-34]。

2.4 挖掘關(guān)聯(lián)分析候選基因

基于中雙11基因組注釋信息,在Bn-A09-p29991443位點下游0.17 Mb和0.18 Mb處分別找到油菜角果長和籽粒質(zhì)量已知基因ARF18和CYP78A9。ARF18編碼1個生長素響應(yīng)因子,CYP78A9編碼1個細(xì)胞色素P450單加氧酶,兩者均通過促進(jìn)油菜角果伸長進(jìn)而影響籽粒質(zhì)量[35-36]。

根據(jù)Darmor基因組的注釋信息,在Bn-A02-p6564861位點上游約0.100 Mb處找到候選基因BnaA02g07770D,其擬南芥同源基因AT5G58230(MSI1)為擬南芥的母性效應(yīng)胚胎滯育基因,在植物胚乳和胚胎發(fā)生過程中起作用,直接影響后期種子的生長發(fā)育,從而導(dǎo)致籽粒質(zhì)量改變[37]。在位點Bn-A06-p19593266上游約0.695 Mb處找到候選基因BnaA06g21460D,其擬南芥同源基因AT5G62000(ARF2)通過限制種皮細(xì)胞增殖影響籽粒大小[38]。在Bn-A01-p15496639位點上游0.635 Mb處找到候選基因BnaA01g20420D,其擬南芥同源基因AT3G48610編碼1個非特異性磷脂酶NPC6,過表達(dá)NPC6能顯著提高亞麻薺的籽粒質(zhì)量、含油量和單株產(chǎn)量[39]。在Bn-A03-p8851142位點上游0.122 Mb處找到候選基因BnaA03g17130D,其擬南芥同源基因AT2G37260(TTG2)編碼的WRKY轉(zhuǎn)錄因子表達(dá)于種皮中,通過增加種皮的細(xì)胞延伸從而促進(jìn)種子的生長發(fā)育[40]。在Bn-A07-p11611255位點上游0.115 Mb處找到候選基因BnaA07g16350D,其擬南芥同源基因AT3G54320(WRI1)編碼1個AP2/ERWEBP類轉(zhuǎn)錄因子,通過增大種子體積和參與貯藏化合物的生物合成以增加種子質(zhì)量[41]。此外,還找到其他16個候選基因,包括擬南芥中已知籽粒質(zhì)量調(diào)控基因CLV3、AHK4、FIS1、GLN1;1、TCTP1、FER、CKX1、EOD1、AHK2等在油菜中的同源拷貝(表3)。

表3 籽粒質(zhì)量候選基因信息

2.5 挖掘WGCNA籽粒質(zhì)量候選基因

利用軟件Cytoscape對purple和magenta模塊中的基因網(wǎng)絡(luò)進(jìn)行可視化表達(dá)(圖5)。在purple模塊中,根據(jù)|KME|≥0.90的標(biāo)準(zhǔn)篩選出340個樞紐基因。樞紐基因BnaA06g00850D(KME=0.97)的擬南芥同源基因AT1G53580(ETHE1)編碼1個硫雙加氧酶,調(diào)節(jié)種子中的硫化物水平,在種子胚胎和胚乳發(fā)育過程中發(fā)揮重要作用[31]。 樞紐基因BnaA01g00990D(KME=0.92)的擬南芥同源基因為籽粒質(zhì)量已知基因AT4G36860(DAR1)(表4),DAR1屬于泛素結(jié)合蛋白家族,通過限制珠被的細(xì)胞增殖從而控制種子大小和籽粒質(zhì)量[42]。樞紐基因BnaC06g10000D(KME=0.91)的擬南芥同源基因AT5G37600(GLN1;1)編碼1種胞質(zhì)谷氨酰胺合成酶,在種子萌發(fā)和灌漿過程中起到?jīng)Q定作用,影響種子胚乳的發(fā)育從而控制籽粒質(zhì)量[43]。此外,purple模塊還包含擬南芥籽粒質(zhì)量已知基因TTN5和SKP1等在油菜中的同源拷貝?；虮磉_(dá)量熱圖(圖6)顯示,這些基因在中雙11中表達(dá)量很低,而在中油821中表達(dá)量較高,特別是在開花后第7 d、10 d和14 d的籽粒和角果皮中,但在成熟的籽粒中表達(dá)量顯著降低,表明它們在籽粒前期發(fā)育中發(fā)揮作用。

表4 Purple和magenta模塊包含的籽粒質(zhì)量已知基因

A:purple模塊內(nèi)基因共表達(dá)網(wǎng)絡(luò);B:magenta模塊內(nèi)基因共表達(dá)網(wǎng)絡(luò)。圖5 Purple和magenta模塊的基因共表達(dá)網(wǎng)絡(luò)Fig.5 Gene co-expression network of purple and magenta modules

a:中油821_角果皮_14 d;b:中油821_角果皮_10 d;c:中油821_角果皮_7 d;d:中油821_種子_45 d;e:中油821_種子_14 d;f:中油821_角果皮_10 d;g:中油821_種子_7 d;h:中雙11_角果皮_14 d;i:中雙11_角果皮_10 d;j:中雙11_角果皮_7 d;k:中雙11_種子_45 d;l:中雙11_種子_14 d;m:中雙11_種子_10 d;n:中雙11_種子_7 d。圖6 Purple和magenta模塊內(nèi)的基因表達(dá)量熱圖Fig.6 Heatmap of gene expression in purple and magenta modules

在magenta模塊中,根據(jù)|KME|≥0.90的標(biāo)準(zhǔn)共篩選到149個模塊樞紐基因,其中樞紐基因BnaC02g44260D(KME=0.96)的擬南芥同源基因為籽粒質(zhì)量已知基因AT5G19400(SMG7)(表4),編碼減數(shù)分裂所必需的蛋白質(zhì)SMG7,SMG7的缺失會引起染色體去密集化和減數(shù)分裂紡錘體的異常重排,從而導(dǎo)致胚死亡,影響籽粒發(fā)育[44],BnaC02g44260D在C09染色體的同源基因BnaC9.SMG7b也已被報道影響油菜每角粒數(shù)、籽粒質(zhì)量等性狀[45]。模塊基因BnaC08g18690D的擬南芥同源基因為籽粒質(zhì)量已知基因AT1G19270(DA1),編碼1個泛素受體,在調(diào)節(jié)胚乳種皮細(xì)胞再復(fù)制中發(fā)揮作用。在甘藍(lán)型油菜中降低DA1的表達(dá)量可以提高種子產(chǎn)量和生物量[46]。模塊基因BnaA06g14090D、BnaA05g14370D和BnaC09g05450D共同對應(yīng)的擬南芥同源基因為AT5G62000(ARF2),通過限制種皮細(xì)胞增殖而影響籽粒大小[38]。此外,magenta模塊還包含擬南芥籽粒質(zhì)量已知基因GLN1;2、RPT2A和BRM等在油菜中的同源拷貝。基因表達(dá)量熱圖(圖6)顯示,這些基因在中油821中表達(dá)量很低,而在中雙11中表達(dá)量較高。在角果皮中,它們的表達(dá)量隨著籽粒發(fā)育進(jìn)程逐漸降低。在籽粒中,這些基因在開花后第7 d和45 d表達(dá)量較高,而在第10 d和14 d表達(dá)量顯著降低,這表明它們在籽粒的發(fā)育早期和末期發(fā)揮作用。

2.6 整合2種分析方法挖掘候選基因

整合GWAS和WGCNA 2種分析方法進(jìn)行分析,結(jié)果表明,在42個顯著關(guān)聯(lián)位點的LD置信區(qū)間內(nèi)找到90個與籽粒質(zhì)量相關(guān)的模塊基因,其中49個屬于magenta模塊,41個屬于purple模塊(表5)。對比它們的擬南芥同源基因功能注釋,發(fā)現(xiàn)部分基因在種子的各項發(fā)育過程中起著決定性作用。在A01染色體位點Bn-A01-p3061223下游0.120 Mb處找到magenta模塊基因BnaA01g06210D,其擬南芥同源基因AT4G30930(NFD1)編碼1個定位在線粒體的核糖體蛋白質(zhì),基因突變后影響雌雄配子的發(fā)育[47]。在A07染色體位點Bn-A02-p305007下游0.170 Mb處找到magenta模塊基因BnaA07g14990D,其擬南芥同源基因AT5G40480(EMB3012)與胚發(fā)育缺陷相關(guān),基因在種子和胚中表達(dá)量較高[48]。在A07染色體位點Bn-A10-p11601681下游0.06 Mb處找到purple模塊基因BnaA07g03030D,其擬南芥同源基因AT2G17800(ROP3)編碼1個GTP結(jié)合蛋白質(zhì),參與生長素激活的Aux/IAA蛋白質(zhì)水解途徑,基因功能突變后導(dǎo)致胚發(fā)育異常[49]。在C03染色體位點Bn-scaff_26320_1-p269450下游0.14 Mb處找到magenta模塊基因BnaC03g46120D,其擬南芥同源基因AT3G01610(CDC48C)參與細(xì)胞分裂循環(huán),與胚發(fā)育缺陷有關(guān)[50]。

表5 顯著關(guān)聯(lián)位點置信區(qū)間內(nèi)的magenta和purple模塊基因

3 討 論

A-D test是一種非參數(shù)統(tǒng)計模型,利用中值替換均值,解決了關(guān)聯(lián)分析中基因型分布不均勻造成的偏差,提高了檢測精度[51]。相較于其他模型,A-D test模型對于效應(yīng)小、數(shù)據(jù)非正態(tài)分布的基因尤其有效。本研究對A-D Test設(shè)置了嚴(yán)格的顯著性閾值(閾值=0.01/總標(biāo)記數(shù)),檢測到19個顯著位點,其中14個位點在至少2個環(huán)境中被重復(fù)檢測到,4個位點在3個環(huán)境中均被重復(fù)檢測到。與另外2個模型的結(jié)果比較,8個位點與GLM重疊,2個位點與MLM重疊。同時,本研究在6個A-D test顯著位點置信區(qū)間內(nèi)找到ARF18、AHK4、FIS1、GLN1;1和ARF2等擬南芥已知籽粒質(zhì)量基因的同源基因。這些結(jié)果證實了A-D test模型的準(zhǔn)確性和可靠性。合并3個模型結(jié)果共得到71個顯著位點,挖掘出ARF2、NPC6、TTG2、WRI1和CLV3等23個擬南芥已知籽粒質(zhì)量基因的同源基因。同時,Bn-A06-p23674542(-lgP=5.97～7.14)等位點在多個環(huán)境和多個模型中被檢測到,但置信區(qū)間內(nèi)尚未找到籽粒質(zhì)量已知基因,后續(xù)可作為研究重點進(jìn)行深入探究。

隨著測序技術(shù)的發(fā)展,利用共表達(dá)網(wǎng)絡(luò)對轉(zhuǎn)錄組進(jìn)行系統(tǒng)研究,是目前比較流行的分析方法。WGCNA能夠聚類表達(dá)模式相似的基因,構(gòu)建與目標(biāo)性狀關(guān)聯(lián)、具有生物學(xué)意義的共表達(dá)模塊,進(jìn)而挖掘共表達(dá)網(wǎng)絡(luò)中的關(guān)鍵基因。在本研究中,我們利用大粒品種中雙11和中小粒品種中油821的28份轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行WGCNA,挖掘到purple和magenta 2個模塊與籽粒質(zhì)量表型顯著相關(guān)。2個模塊基因顯著富集條目包括苯丙氨酸氨基轉(zhuǎn)移酶活性、硫雙加氧酶活性、肉桂酸生物合成和焦磷酸酶活性等,這些酶活性在種子生長發(fā)育過程中都起到重要作用。2個模塊的樞紐基因如BnaA06g00850D、BnaA01g00990D、BnaA07g36220D和BnaC06g10000D等擬南芥同源基因ETHE1、DAR1、TTN5和GLN1;1等已被報道參與籽粒質(zhì)量調(diào)控,magenta模塊樞紐基因BnaC02g44260D是油菜已克隆每角粒數(shù)、籽粒質(zhì)量調(diào)控基因BnaC9.SMG7b在C02染色體的同源拷貝。這些結(jié)果證實了通過WGCNA挖掘油菜籽粒質(zhì)量性狀調(diào)控基因的可行性和可靠性。

GWAS分析的難點之一是在顯著位點置信區(qū)間內(nèi)高效鑒定候選基因,而WGCNA能夠特異地篩選出與目標(biāo)性狀關(guān)聯(lián)的差異表達(dá)基因,2種分析方法具有很好的互補(bǔ)性。Farber[52]首次結(jié)合GWAS和WGCNA挖掘調(diào)控骨密度的相關(guān)基因,對關(guān)聯(lián)候選基因構(gòu)建共表達(dá)網(wǎng)絡(luò),挖掘到3個與性狀顯著相關(guān)的基因模塊。在油菜的研究中,越來越多的研究者也采取相同策略,如王艷花等[53]結(jié)合GWAS和WGCNA 2種方法篩選出調(diào)控油菜生物產(chǎn)量的候選基因BnA07g19320D、BnC03g33610D和BnC03g73810D等,鮮小華等[20]結(jié)合GWAS和WGCNA 2種方法挖掘出甘藍(lán)型油菜黃籽候選基因BnATCAD4、BnF3H和BnANS等。在本研究中,我們在42個顯著位點的置信區(qū)間內(nèi)找到90個屬于2個目標(biāo)模塊(purple模塊和magenta模塊)的基因。盡管這些基因尚未被報道直接參與籽粒質(zhì)量調(diào)控,但它們在大粒品種中雙11和中小粒品種中油821的籽粒中表達(dá)豐度不同,部分基因如BnaA01g06210D、BnaA07g14990D、BnaA07g03030D等直接或間接調(diào)節(jié)種子、胚胎等相關(guān)發(fā)育進(jìn)程,后續(xù)工作可以深入挖掘這些基因在籽粒質(zhì)量調(diào)控中的具體功能。

4 結(jié) 論

本研究利用A-D Test模型對496份油菜的籽粒質(zhì)量進(jìn)行GWAS分析,共檢測到19個顯著位點。整合前期MLM和GLM模型結(jié)果后共得到71個位點,聯(lián)合解釋50.1%的表型變異。在21個關(guān)聯(lián)位點置信區(qū)間內(nèi)找到SMG7、ARF2、NPC6、TTG2、WRI1等油菜、擬南芥中已被報道的籽粒質(zhì)量基因的同源基因。利用大粒品種中雙11和中小粒品種中油821的28份轉(zhuǎn)錄組數(shù)據(jù),構(gòu)建了13個共表達(dá)模塊,其中magenta和purple模塊與籽粒質(zhì)量表型相關(guān)性最高。在2個目標(biāo)模塊中挖掘出BnaA06g00850D、BnaA01g00990D、BnaC06g10000D、BnaC02g44260D和BnaC08g18690D等籽粒質(zhì)量候選基因。整合GWAS和WGCNA的分析結(jié)果,在42個顯著位點的置信區(qū)間內(nèi)找到90個屬于2個籽粒質(zhì)量相關(guān)的模塊基因,其中BnaA01g06210D、BnaA07g14990D和BnaA07g03030D等參與調(diào)控種子、胚胎等相關(guān)發(fā)育進(jìn)程;同時,A-D Test挖掘的多個大效應(yīng)位點如Bn-A06-p23674542等置信區(qū)間內(nèi)尚未找到與籽粒質(zhì)量相關(guān)的已知基因,值得后續(xù)深入挖掘。