基于轉錄組測序的關聯分析定位裂葉樺葉形調控基因

顧宸瑞,袁啟航,姜 靜,穆懷志,劉桂豐*

(1. 林木遺傳育種全國重點實驗室(東北林業大學),黑龍江 哈爾濱 150040;2. 北華大學林學院, 吉林 吉林 132013)

作為歐洲白樺(Betulapendula)的變種,裂葉樺(B.pendula‘Dalecarlica’)葉片開裂成掌狀,頂部葉尖細長,一級齒減少,葉脈在葉背面明顯隆起,內部機械輸導組織非常發達,且橫切面積顯著增大,葉背絨毛發達,在園藝觀賞、育種工作及植物發育生物學研究中具有重要的應用價值,但其葉片形態變異機制及調控基因始終未被揭示[1-2]。

雙子葉植物的葉片性狀呈現非常豐富的多樣性,一般從葉形指數、葉緣缺刻程度和葉緣鋸齒3個方面對葉片的形狀進行衡量,而這3個指標分別受葉側脈延伸速率和葉肉蹼化程度兩個因素的控制,葉側脈延伸速率大則葉片較寬,葉肉的蹼化程度在整個葉緣較低則出現鋸齒,蹼化過程在某幾個葉脈間的位置受抑制則出現掌狀裂[3]。

植物的葉形可能是數量性狀或質量性狀,其受環境影響嚴重,且某些植物存在異形葉性,決定葉形的成因復雜,這給定位控制葉形性狀的主效應基因帶來一定難度[4]。目前對雙子葉植物葉片形狀的多樣性及其遺傳規律的研究主要集中于棉花(Gossypiumhirsutum)、羽衣甘藍(Brassicaoleraceavar.acephala)、白菜(Brassicarapevar.glabra)等葉面積能夠直接或間接影響產量的作物上。其中,影響棉花葉片形狀的基因已被定位于其第15號染色體[5],GhLMI1-like和GhKNOX1基因對葉形的控制已被驗證[6];羽衣甘藍的葉緣缺刻性狀被認為是數量性狀,但起決定性的主效應基因共有2對[7-9];白菜不結球裂葉突變體的裂葉性狀被發現是受一對顯性基因控制的質量性狀,并已發現與之關聯的遺傳標記[10]。

關聯分析以連鎖不平衡為基礎,利用基因組上的分子標記,定位特定性狀對應的基因座。隨著高通量測序技術的發展,全基因組關聯分析(GWAS)已成為研究群體數量性狀變異的常用手段[11-14]。GWAS策略利用規模較大的野生群體和全基因組高密度SNP標記定位目標基因座,可節省群體構建成本和時間。但這種關聯分析需要對目標性狀進行精確的測量和量化,而葉形不規則程度(復葉,裂片、葉緣鋸齒等特征)缺乏合適的量化標準,對這些數量性狀的關聯分析幾乎沒有報道;另一方面,GWAS使用DNA重測序數據,不能同時檢測基因的變異和表達量差異,而植物葉片的形態建成往往受基因表達量影響[15]。轉錄組測序可一次性獲取外顯子區的SNP標記和基因表達量信息,但受限于無法檢測基因間區SNP標記的缺點,基于轉錄組測序的關聯分析未被廣泛采用。

近年來,在裂葉樺與歐洲白樺的雜交實驗中,裂葉樺葉形始終被作為質量性狀進行統計[2],但孟德爾分離規律不能解釋子代群體中出現的中間過渡型個體和超親遺傳個體。本研究提出了一種量化葉形不規則程度的指標,計量了24株不同葉形白樺的葉不規則程度,并與其RNA-Seq數據進行關聯分析,旨在定位控制白樺葉形的主效應基因座,預測候選基因,為后續的驗證和育種工作奠定基礎。

1 材料與方法

1.1 供試材料

試驗材料來自裂葉樺(中裂葉)×歐洲白樺的一個半同胞家系子代,定植于東北林業大學白樺強化種子園。該家系有卵形葉、深裂葉、中裂葉和淺裂葉4種葉形的后代(圖1),每種選取6株,共24株,用于后續研究。

圖1 實驗材料不同類型葉片示例Fig. 1 Leaf shape examples for experimental material

1.2 試驗方法

1.2.1 葉不規則度的測量和比較

深裂葉類型每株摘取3片成熟功能葉,其他葉形每株隨機摘取5片發育成熟的功能葉,共108片葉。使用Epson Expression 10000XL掃描儀(日本)獲取葉片圖像,所有圖片以300 dpi進行掃描。

1.2.2 影響裂葉樺葉形基因的定位

24個植株頂端分生組織RNA-Seq數據來自本團隊的前期工作[16]。使用RNA-STAR軟件[17]將測序數據比對到歐洲白樺參考基因組(CoGe:35080)[18],使用sambamba軟件[19]去除重復讀段,使用GATK4分析SNP/InDel,并將24個植株的SNP/InDel信息匯總。使用beagle軟件[20-21]進行數據填補。使用plink軟件[22]進行過濾,標準為:MAF≥5%,GENO≤20%,MIND≤20%,PHWE≥0.001(MAF表示第2等位基因頻率,GENO表示位點基因型缺失比例,MIND表示樣本基因型缺失比例,PHWE表示哈迪溫伯格平衡顯著性)。使用Tassel 5.0軟件[23]進行基于一般線性模型(GLM)的關聯分析。以P=0.01/n計算Bonferroni 閾值(n為參與關聯分析的SNP/InDel個數)篩選目標SNP。

1.2.3 候選基因及注釋

根據參考基因組注釋文件查找目標SNP覆蓋的基因。使用htseq-count軟件[24]依據RNA-Seq數據計算全基因組基因表達量,使用DESeq2軟件[25]計算基因表達差異顯著性(將來自6株正常葉個體的數據視為野生型的6次生物學重復,將來自18株裂葉個體的數據視為裂葉樺的18次生物學重復),并篩選出位于關聯區的差異基因。使用COG、Gene Ontology、KEGG、KOG、Pfam、Swissprot、TrEMBL、NCBI-nr等在線數據庫對候選基因進行注釋。

2 結果與分析

2.1 實驗材料表型數據

根據裂葉樺×歐洲白樺(B.pendula)半同胞家系24株子代共108片成熟葉的圖像計算葉不規則程度,并進行方差分析(表1)和多重比較(表2)。

表1 葉形不規則度方差分析結果

表2 個體間葉形不規則程度多重比較結果

由表1可見,葉不規則程度在不同個體之間差異極顯著(P<0.01)。由表2可見,葉形不規則程度變幅較大,不規則指數最高可達0.352 2,最低僅0.032 8,平均0.182 4。樣本變異系數為60.51%,滿足進行關聯分析的條件。多重比較結果表明,該量化指標在不同葉形個體間差異顯著,其數值能夠代表不同個體的葉形。

2.2 關聯分析及候選基因的預測與注釋

RNA-Seq數據系來自Illumina HiSeq 2500平臺的雙端100 bP第2代高通量測序結果,共331 Gb,平均每個樣品13.8 Gb。過濾后,Q30比例均在90.62%以上,與參考基因組比對率均在85.36%以上。

由RNA-Seq數據共獲得SNP/InDel 3 303 530個,經過補齊和篩選,保留了高質量多態SNP/InDel 1 367 359個。采用GLM模型進行關聯分析,發現超過Bonferroni 閾值(-lgP≥ 8.135 883)的候選SNP/InDel共66個(圖2 a、2 b)。由于實驗材料來自同一半同胞家系,將零星出現于3、6、11號染色體上的5個候選SNP視為異常值剔除。

a. 曼哈頓圖Manhattan plots; b.QQ圖 Quantile-quantile plots;c.7號染色體上目標SNP在個體中的基因型 genotypes of aim SNPs on chromosome7。圖2 葉片不規則程度關聯分析結果Fig. 2 Association analysis result of leaf irregularity degree

保留的61個目標SNP/InDel分布于一個9.19 MbP的關聯區中(BpChr07:7466344—16651477)。該區域共包含400個基因,其中17個基因的編碼區攜帶上述目標SNP/InDel。位于基因編碼區的SNP基因型均不能與表型完全關聯,但這些SNP普遍在正常葉個體與裂葉個體之間基因型存在差異,除深裂葉個體C-3外,目標SNP在裂葉樺中基因型一致(圖2 c)。

24個樣品的全基因組基因表達量及差異顯著性由DESeq2分析,在關聯區中,表達量在裂葉個體與正常葉個體之間存在極顯著差異(P<0.01)的基因共25個,其中6個在裂葉樺中上調表達,19個在裂葉樺中下調表達(圖3 a)。

被目標SNP影響編碼的基因和關聯區差異表達基因共計40個,其中2個同時具有SNP和表達量差異(圖3 b)。使用公共數據庫對候選基因進行了功能注釋,在表3中,對注釋結果進行了概括性的描述。

表3 候選基因的變異類型及注釋信息

由表3可見,候選基因中包括轉錄因子3個、作為細胞組成成分的蛋白3個、參與蛋白質翻譯加工的蛋白6個、參與信號傳導的蛋白激酶9個、參與生化反應的酶5個、參與物質運輸的蛋白3個、對RNA進行轉錄后修飾的蛋白1個,DNA內/外切酶1個,以及調節氣孔的蛋白1個。另外,8個基因無法在現有的在線數據庫中查詢到,其類型與功能完全未知。

3 討 論

3.1 Dli(4πS/C2)可用作量化植物葉形不規則程度的指標

植物葉片邊緣的形態通常以“鋸齒”“裂葉”“深裂”“復葉”等非量化指標進行描述[15],而缺少將葉片形態的不規則程度進行量化的指標,因此,植物的葉形變異通常被作為質量性狀處理,或被解釋為異形葉性。本研究所使用的半同胞家系,個體出現從卵形葉到雞爪狀深裂葉的分布,為了進行關聯分析,需要對葉緣開裂程度進行量化。本研究以Dli作為衡量葉片不規則程度的指標,該指標為平面圖形面積與周長相等的圓的面積之比。

值得注意的是,在計算葉片圖像周長時,會受到尺縮效應的影響,即圖片精細度不足以展示葉緣鋸齒時,會使測得的周長偏小;而葉緣絨毛、破損等細節被過度放大時,會使測得的周長偏大。因此對每一批次的比較,必須先通過預實驗選定適用于該物種的像素密度,而后在相同的像素密度下掃描樣本圖像。

3.2 基于轉錄組數據進行關聯分析的優勢與不足

隨著高通量測序技術的發展,以GWAS為代表的利用高密度SNP為標記的關聯分析方法被廣泛應用于數量性狀基因座的定位。通常,GWAS所需的樣本來自大規模的野生群體,數量在100以上,樣本的數量性狀表型值呈連續的正態分布[26]。

本研究中,對24株來自同一半同胞家系的材料進行轉錄組測序,并通過一般線性模型(GLM)進行了關聯分析。由于實驗材料來自一個半同胞家系,且樣本數量較少,這使得實際P值過低,高于閾值的SNP過多,觀測-lgP值總體偏高,不服從正態分布,因而不能簡單地套用GWAS的分析策略。本研究表明,在半同胞家系中,相鄰的SNP標記體現出連鎖效應。筆者效仿BSR分析的策略[27],引入關聯區的概念,一方面將關聯區外零星出現的高于閾值的SNP視為異常值剔除,另一方面對關聯區中所有發生了可能與表型有關的變異(包括關聯性SNP/InDel與差異表達)的基因進行注釋。結果表明,調控裂葉樺葉形的基因坐落于一個9.18MbP的關聯區,其中40個基因出現了可與表型相關聯的編碼變異或表達量差異。

為了避免第一類錯誤,筆者匯總了全部40個在不同葉形個體之間存在差別的基因,統一進行了基因功能注釋。但本研究未檢測到與表型完全關聯的SNP,這可能是由于裂葉樺的葉形是多個基因突變的綜合調控結果,或由某基因的啟動子變異導致。能夠一次性獲得序列變異和表達量差異是RNA-Seq相較于DNA-Seq的優勢,但RNA-Seq不能檢測基因啟動子的序列,因此對差異表達的基因,應在后續研究中關注其啟動子是否存在變異。

3.3 關聯分析結果為基因精細定位奠定了基礎

候選基因中共包含3個轉錄因子,其中2個隸屬于MYB家族(Bpev01.c1499.g0007、Bpev01.c1097.g0004),1個隸屬于WOX家族(Bpev01.c1525.g0003)。在擬南芥(Arabidopsisthaliana)中,AS1作為MYB結構域蛋白,可通過與KNOX等轉錄因子拮抗或互作調控葉片近/遠軸的形態建成從而影響葉形[28];在煙草(Nicotianasylvestris)中,WUS作為WOX家族成員,參與分生組織的維持,是葉原基形成和葉脈建成所必需的[29]。候選基因中的一個PIN家族蛋白(Bpev01.c0147.g0002)在白樺(B.platyphylla)[30]中的同源基因BpPIN1已被克隆,BpPIN1的表達量與白樺葉片發育過程緊密相關,并在白樺與裂葉樺之間存在顯著表達模式差異[31]。在擬南芥中,PIN1作為介導生長素極性運輸的外泌蛋白,直接制造葉原基內部特定位置的生長素濃度峰點,繼而誘導這些位置發育成葉脈[32]。裂葉樺裂葉性狀的成因可能與上述4個候選基因有關,但不能排除其他候選基因參與調控葉形的可能,尤其是8個未能在現有數據庫中查詢到的候選基因。這8個未知基因可能是白樺所特有的,暫無法對其功能進行預測。

本研究所用材料來自一個半同胞家系,由于此類近緣群體中相鄰SNP的連鎖較強,導致定位結構所得關聯區較大。使用近緣群體進行關聯分析的同類研究案例表明,此類關聯分析是可行的,其結合精細定位工作成功確定了目的基因[33]。精細定位可縮小關聯區范圍,大量地排除因近緣群體SNP連鎖及隨機陽性信號帶來的第二類錯誤,排除與葉形無關的候選基因。在后續研究中,可使用分子生物學手段,借助基因組上已發掘的SSR標記[18, 30]和本研究中檢測到的SNP標記,在本研究基礎上開展精細定位工作,將關聯區進一步縮小。