劉嘉智，劉長寧

（1.中國科學(xué)院西雙版納熱帶植物園，云南勐臘 666303；2.中國科學(xué)院大學(xué)，北京 100049）

橡膠籽是天然橡膠種植生產(chǎn)中的重要副產(chǎn)品，橡膠籽的含油量為35%～50%，年產(chǎn)量可達(dá)2 060 kg/hm2［1］，與其 他 非 食用 植 物 油相 比，橡膠籽油燃燒的CO2和NO2排放率更低，是理想的生物柴油來源［2-3］。然而，目前對橡膠樹的利用主要集中于獲取天然橡膠，忽略了橡膠籽作為生物柴油的巨大潛力，與之相關(guān)的研究也并不充足。

通常，具有同源序列的蛋白編碼基因具有相似的功能，因而基于同源分析可以在近緣物種中預(yù)測相關(guān)功能基因。本研究基于共線性分析和同源分析，利用目前研究中較為深入的蓖麻全基因組關(guān)聯(lián)分析數(shù)據(jù)，將蓖麻野生型品種中種子農(nóng)藝性狀相關(guān)的基因映射到橡膠樹reyan7-33-97 品種基因組上，使用生物信息學(xué)手段預(yù)測橡膠樹中種子性狀相關(guān)基因并進(jìn)行基因功能分析，為橡膠籽的合理開發(fā)利用提供理論參考。

1 材料和方法

1.1 數(shù)據(jù)來源

蓖麻野生型品種的基因組文件以及基因組注釋文件下載于http://oilplants.iflora.cn/Download/castor_download.html，橡膠樹reyan7-33-97 品種的基因組文件以及基因組注釋文件從https://www.ncbi.nlm.nih.gov/data-hub/genome/GCF_001654055.1/下載。從上述文件分別提取基因ID、序列信息、基因組位置信息和功能描述等內(nèi)容。從NCBI SRA 數(shù)據(jù)庫中收集具有文獻(xiàn)來源且為雙端測序的蓖麻和橡膠樹根、莖、種子3 種組織的轉(zhuǎn)錄組數(shù)據(jù)，樣本編號分別為SRR3136156、SRR6127585、SRR3136168、SRR17218208、SRR17218209、SRR17218216。

1.2 研究方法

1.2.1 橡膠樹種子性狀相關(guān)基因預(yù)測

根據(jù)文獻(xiàn)調(diào)研［4］，在蓖麻野生型品種中收集基于全基因組關(guān)聯(lián)分析得到的種子特征性狀相關(guān)的SNP 位點(diǎn)。根據(jù)基因組注釋信息，在上述SNP 位點(diǎn)前后50 kb 范圍內(nèi)提取基因，并認(rèn)為這些基因可能與種子特征有關(guān)。

橡膠樹與蓖麻同屬大戟科，具有較近的親緣關(guān)系，因此，首先判斷蓖麻種子性狀相關(guān)基因所在區(qū)塊在橡膠樹中是否存在共線區(qū)塊。由于橡膠樹reyan7-33-97 品種基因組只組裝到scaffold 級別，使用BLAST 進(jìn)行共線性分析后認(rèn)為存在3 個(gè)及以上基因，且基因之間相互間隔不超過50 kb的片段具有共線性。然后，使用Orthofinder 軟件［5］對這兩個(gè)物種進(jìn)行同源組鑒定，得到的結(jié)果再根據(jù)雙向BLAST 最優(yōu)策略進(jìn)一步確定同源基因，并根據(jù)基因結(jié)構(gòu)驗(yàn)證預(yù)測的同源基因是否可靠，最終將收集到的蓖麻野生型種子性狀相關(guān)基因比對到橡膠樹基因組上，來挖掘橡膠樹中與種子特征性狀相關(guān)的基因。

1.2.2 表達(dá)譜分析和差異表達(dá)分析

下載蓖麻和橡膠樹根、莖、種子3 種組織對應(yīng)的轉(zhuǎn)錄組原始測序文件，使用Trim Galore 進(jìn)行質(zhì)控檢測去除低質(zhì)量堿基和接頭。利用HISAT2 軟件［6］對基因組序列構(gòu)建索引后進(jìn)行序列比對，使用Stringtie 工具［7］進(jìn)行定量分析來獲取基因水平上的Raw counts、TPM 值以及FPKM 值。將TPM 值進(jìn)行歸一化處理，使用R 包pheatmap 繪制表達(dá)量熱圖，并對上述兩個(gè)物種得到的表達(dá)量結(jié)果進(jìn)行比較。使用GFOLD 1.1.4 軟件［8］對橡膠樹基因進(jìn)行差異表達(dá)分析，進(jìn)一步挖掘橡膠樹中與種子性狀相關(guān)基因在種子與根、莖之間的表達(dá)差異。

1.2.3 miRNA 靶位點(diǎn)預(yù)測

根據(jù)文獻(xiàn)調(diào)研［9］以及pmiren 數(shù)據(jù)庫［10］中獲取橡膠樹的成熟miRNA 序列，使用gffread 軟件提取轉(zhuǎn)錄本序列，使用psRNATarget 在線工具［11］默認(rèn)參數(shù)進(jìn)行miRNA 靶基因預(yù)測，從預(yù)測結(jié)果中篩選期望值大于3.5 的相互作用的miRNA 和靶基因，使用Cytoscape 軟 件［12］將miRNA 與 靶 基因 間的 網(wǎng)絡(luò)關(guān)系進(jìn)行可視化展示。

1.2.4 蛋白互作網(wǎng)絡(luò)預(yù)測

基于同源組劃分，以擬南芥蛋白互作網(wǎng)絡(luò)作為參考映射到橡膠樹中，獲得橡膠樹reyan7-33-97 的蛋白互作網(wǎng)絡(luò)。從String 數(shù)據(jù)庫［13］中下載擬南芥蛋白互作物理子網(wǎng)，使用Orthofinder 軟件預(yù)測擬南芥與橡膠樹reyan7-33-97 之間的蛋白同源群組。使用Cytoscape 軟件對橡膠樹種子性狀相關(guān)候選基因進(jìn)行網(wǎng)絡(luò)可視化。

2 結(jié)果與分析

2.1 橡膠樹種子性狀相關(guān)基因預(yù)測

在野生型蓖麻中共發(fā)現(xiàn)18 個(gè)與種子特征性狀相關(guān)的SNP 位點(diǎn)，這些SNP 位點(diǎn)分別與種子面積（SA）、種子長度（SL）、種子厚度（ST）、種子寬度（SW）以及單粒重（SSW）等5 個(gè)特征形狀相關(guān)。根據(jù)蓖麻基因組注釋文件，在上述SNP 位點(diǎn)上下50 kb 范圍內(nèi)共有67 個(gè)基因（圖1），分別位于蓖麻的1、2、3、5、6 以及9 號染色體上（表1）。其中Rc01G001599、Rc01G001600、Rc01G001601、Rc01 G001602、Rc01G001603、Rc01G001604、Rc01G001605、Rc01G001606、Rc01G001607、Rc01G001608、Rc01G 001609 對上述5 個(gè)性狀都產(chǎn)生影響。

表1 蓖麻野生型種子特征相關(guān)基因

圖1 蓖麻野生型種子性狀相關(guān)基因統(tǒng)計(jì)

基于共線性分析，在蓖麻野生型與橡膠樹reyan7-33-97 中發(fā)現(xiàn)4 個(gè)共線區(qū)塊，分別位于蓖麻的2、5、6 以及9 號染色體上，對應(yīng)橡膠樹的NW_018746954.1、NW_018745856.1、NW_018746236.1 以及NW_018746078.1 這4 條scaffold，共 覆 蓋19 個(gè)基因（表2）。根據(jù)Orthofinder 軟件，在蓖麻野生型與橡膠樹reyan7-33-97 中共鑒定到16 926 個(gè)同源組，根據(jù)雙向BLAST 最優(yōu)確定二者之間序列相似程度，上述67 個(gè)蓖麻基因最終在橡膠樹中比對得到49 個(gè)基因。對上述橡膠樹基因以及同源蓖麻基因根據(jù)基因結(jié)構(gòu)作圖進(jìn)行驗(yàn)證（圖2），發(fā)現(xiàn)基因結(jié)構(gòu)與最終得到的基因序列預(yù)期結(jié)果相符。

表2 蓖麻野生型與橡膠樹reyan7-33-97共線區(qū)塊基因統(tǒng)計(jì)

圖2 蓖麻野生型與橡膠樹reyan7-33-97同源基因結(jié)構(gòu)圖

2.2 種子性狀相關(guān)表達(dá)譜和差異表達(dá)分析

為進(jìn)一步探究蓖麻和橡膠樹種子性狀相關(guān)基因的功能，進(jìn)行了組織表達(dá)譜的分析，圖3 是蓖麻和橡膠樹種子性狀相關(guān)基因的組織表達(dá)譜。從中發(fā)現(xiàn)，有些同源基因在3 個(gè)組織中的表達(dá)量類似，但有些差別很大。對蓖麻和橡膠對應(yīng)的49 個(gè)基因在種子中表達(dá)結(jié)果進(jìn)行統(tǒng)計(jì)，發(fā)現(xiàn)在兩個(gè)物種中有5 個(gè)基因均為正表達(dá)，25 個(gè)基因均為負(fù)表達(dá)，19 個(gè)基因表達(dá)相反，這可能是由于物種在分化過程中部分基因功能改變導(dǎo)致的。

圖3 蓖麻和橡膠樹種子性狀相關(guān)基因組織表達(dá)譜

根據(jù)差異表達(dá)分析，選出16 個(gè)在種子中表達(dá)量遠(yuǎn)高于根莖的基因、17 個(gè)在種子中表達(dá)量遠(yuǎn)低于根莖的基因以及5 個(gè)在三者之中都沒有差異表達(dá)的基因，其余基因在3 個(gè)組織中無交集（表3）。

表3 種子與根莖差異表達(dá)基因

2.3 miRNA 靶位點(diǎn)預(yù)測

基因的miRNA 靶位點(diǎn)預(yù)測也是研究基因調(diào)控的重要途徑。本研究共收集到橡膠樹的217 條成熟miRNA 序列，對橡膠樹種子性狀相關(guān)基因進(jìn)行miRNA 靶位點(diǎn)預(yù)測，結(jié)果發(fā)現(xiàn)LOC110671483、LOC110673189、 LOC110645416、 LOC110651200、LOC110672934、 LOC110641727、 LOC110653311、LOC110659520 受到不同miRNA 調(diào)控（圖4）。由于miRNA 在物種進(jìn)化中相當(dāng)保守，大部分miRNA 在特定的組織和發(fā)育階段表達(dá)，我們預(yù)測得到的miRNA 許多功能都與種子性狀相關(guān)。例如之前研究發(fā)現(xiàn)miR156 參與調(diào)控?cái)M南芥種子成熟過程［14］、miR408 在水稻種子發(fā)育中具有重要作用［15］等，因此推測上述8 個(gè)基因可能受到miRNA 調(diào)控影響種子發(fā)育，從而影響了種子特征性狀。

圖4 橡膠樹種子性狀相關(guān)基因miRNA靶位點(diǎn)預(yù)測

2.4 蛋白互作網(wǎng)絡(luò)預(yù)測

基于同源組的劃分，橡膠樹reyan7-33-97 品種中共有22 815 個(gè)基因映射到19 762 個(gè)擬南芥基因。根據(jù)預(yù)測得到的橡膠樹蛋白互作網(wǎng)絡(luò)，我們將橡膠樹reyan7-33-97 品種中與種子性狀相關(guān)的49 個(gè)基因放入蛋白互作網(wǎng)絡(luò)中，發(fā)現(xiàn)除LOC110654127 和LOC110668244 外，其 它 基 因 之 間不存在邊。這可能是影響種子性狀的橡膠樹候選基因與其他基因也存在相互作用導(dǎo)致的。為此，提取這49 個(gè)基因以及最近的相鄰基因重新進(jìn)行網(wǎng)絡(luò)構(gòu)建，結(jié)果發(fā)現(xiàn)其中16 個(gè)橡膠樹種子性狀候選基因通過一級相互連接構(gòu)成網(wǎng)絡(luò)（圖5），因此認(rèn)為這16 個(gè)基因存在間接蛋白相互作用從而影 響 種 子 性 狀，LOC110654127 和LOC110668244 同時(shí)存在直接蛋白相互作用影響種子性狀，而其它33 個(gè)基因可能單獨(dú)影響種子性狀或者存在不同于擬南芥的蛋白互作成員。

圖5 橡膠樹種子性狀候選基因蛋白互作網(wǎng)絡(luò)

3 展望

本研究使用全基因組關(guān)聯(lián)分析、基因組共線性分析和基因序列同源分析等多種技術(shù)手段，預(yù)測了49 個(gè)與橡膠樹種子性狀相關(guān)的基因并進(jìn)行了多種組學(xué)分析，研究結(jié)果為進(jìn)一步開發(fā)利用橡膠籽提供了一定的工作基礎(chǔ)。然而，橡膠籽種子性狀基因與種子產(chǎn)油之間的關(guān)系是一個(gè)復(fù)雜的過程，植物激素的生物合成和信號傳遞對種子大小和形態(tài)的控制至關(guān)重要，例如赤霉素可以促進(jìn)種子的擴(kuò)大和發(fā)育［16］，而乙烯則可以抑制種子的生長［17］。因此，后續(xù)研究可深入挖掘橡膠籽種子性狀基因和植物激素生物合成和信號傳遞通路之間的關(guān)系，通過對橡膠樹中這些激素合成和信號傳遞通路進(jìn)行深入分析，揭示橡膠籽種子性狀基因在種子性狀形成過程中的作用機(jī)制，以便更好地理解種子大小和形態(tài)的調(diào)控機(jī)制，為橡膠樹的育種和生產(chǎn)提供指導(dǎo)。