桑樹(shù)基因組研究團(tuán)隊(duì)

（家蠶基因組生物學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，重慶 400716）

桑樹(shù)為落葉喬木，對(duì)家蠶而言是一種重要經(jīng)濟(jì)飼料作物。桑葉飼養(yǎng)家蠶用于生產(chǎn)貴重蠶絲的過(guò)程至少起始于5 000年以前［1］，通過(guò)“絲綢之路”塑造了世界歷史。

桑科包含了37個(gè)屬，約1 100個(gè)種，其中包括大家熟知的桑樹(shù)、面包樹(shù)、無(wú)花果樹(shù)、榕樹(shù)、見(jiàn)血封喉樹(shù)等［2］。桑樹(shù)屬于桑屬，桑屬包含了10到13個(gè)公認(rèn)的種，超過(guò)1 000個(gè)栽培品種［3］。這些栽培品種廣泛的種植于歐亞大陸、非洲以及美國(guó)。在中國(guó)和印度，分別有約626 000、280 000hm2的土地用來(lái)栽桑養(yǎng)蠶［4］。由于美味的桑椹、造紙用的桑皮及其在東方傳統(tǒng)醫(yī)學(xué)方面的多種用途等特性，桑樹(shù)的栽植對(duì)農(nóng)民具有很強(qiáng)的吸引力［5－6］。

家蠶是一個(gè)專食桑葉的昆蟲(chóng)，為鱗翅目的模式昆蟲(chóng)。目前已知的大多數(shù)鱗翅目昆蟲(chóng)都是農(nóng)業(yè)及林業(yè)上的害蟲(chóng)，因此對(duì)經(jīng)濟(jì)具有重要的影響。促攝食物質(zhì)對(duì)理解植物－昆蟲(chóng)相互之間關(guān)系具有關(guān)鍵的作用，家蠶成功馴化使得深入研究這些促攝食物質(zhì)成為可能，這對(duì)理解植物與昆蟲(chóng)之間關(guān)系是非常關(guān)鍵的。家蠶基因組精細(xì)圖已于2008年測(cè)序完成［7－8］。然而，桑屬中物種幾乎沒(méi)有任何基因組信息。桑樹(shù)基因組測(cè)序完成不僅能夠促進(jìn)桑樹(shù)改良，而且通過(guò)對(duì)桑樹(shù)與家蠶基因組的分析將進(jìn)一步加深我們對(duì)植物－植食性昆蟲(chóng)之間適應(yīng)性進(jìn)化的理解。

本文我們報(bào)道了桑樹(shù)（川桑）基因組測(cè)序草圖。桑樹(shù)基因組估算為357MB，包含了7對(duì)染色體，通過(guò)Illumina技術(shù)測(cè)序覆蓋深度達(dá)到236倍?；诮M裝完成的330MB基因組，鑒定到128MB的重復(fù)序列以及29 338個(gè)蛋白編碼基因。比較基因組分析表明相較于其他已測(cè)序的薔薇目植物而言，桑樹(shù)進(jìn)化速度更快。對(duì)桑樹(shù)抗性基因的鑒定及分析將加速桑樹(shù)品種改良。另外，我們?cè)诩倚Q兩個(gè)組織中發(fā)現(xiàn)了桑樹(shù)的miRNA，這個(gè)結(jié)果從分子水平上暗示了植物－植食性昆蟲(chóng)之間可能的相互關(guān)系。

1 結(jié) 果

1.1 基因組測(cè)序和組裝

川桑體細(xì)胞包含7對(duì)不同的染色體（圖1），我們通過(guò)全基因組鳥(niǎo)槍法策略對(duì)川桑進(jìn)行測(cè)序。合計(jì)78.34billion高質(zhì)量的堿基（覆蓋深度236倍）組裝成330.79MB的桑樹(shù)基因組。Scaffold N50長(zhǎng)度為390 115bp，重疊群N50長(zhǎng)度為34 476bp（表1、附表1和附表2）。組裝完成的基因組中包含16 281kb（4.9%）的gaps以及314 510kb（95.1%）的無(wú)缺口連續(xù)序列。我們從500bp插入片段文庫(kù)中選擇了10.46Gb的高質(zhì)量reads用來(lái)計(jì)算K－mer深度（此處選擇17bp）的分布。總共得到8 577 674 309的17－mer，據(jù)此估計(jì)川桑的基因組為357.4MB（附方法、附圖1和附表3）。681個(gè)scaffolds覆蓋了超過(guò)80%的組裝序列，其中最長(zhǎng)的scaffold為3 477 367bp，93.96%的堿基覆蓋深度超過(guò)20reads（附圖2），在隨機(jī)選擇的10 000個(gè)表達(dá)序列標(biāo)簽中，97%的表達(dá)序列標(biāo)簽長(zhǎng)度的90%被一個(gè)scaffold覆蓋（附圖4）。桑樹(shù)基因組GC含量為35.02%，與其他雙子葉植物GC含量類(lèi)似（附方法和附圖3）。

圖1 川桑（M.notabilis）染色體細(xì)胞學(xué)分析

表1 桑樹(shù)基因組測(cè)序及組裝數(shù)據(jù)統(tǒng)計(jì)

1.2 重復(fù)序列

在無(wú)缺口的桑樹(shù)基因組中結(jié)合從頭預(yù)測(cè)以及基于Repbase文庫(kù)同源搜索的方法共得到了127.98MB的重復(fù)序列（附表5）。由于從頭測(cè)序技術(shù)在處理重復(fù)序列方面固有的限制，桑樹(shù)基因組中轉(zhuǎn)座元件的含量有可能被低估。在桑樹(shù)基因組中根據(jù)平均覆蓋深度、總reads所覆蓋到重復(fù)序列（約127.7MB）以及非重復(fù)序列區(qū)域（約166.0MB），排除未測(cè)出的堿基（N）后，我們估計(jì)在未組裝完的序列中約有18.48MB的重復(fù)序列。因此，在桑樹(shù)基因組中存在約有47%的重復(fù)序列。這個(gè)比例與蘋(píng)果中重復(fù)序列的比例（42%）比較接近，略高于楊樹(shù)（35%）。超過(guò)50%的桑樹(shù)重復(fù)序列可以劃分到已知的類(lèi)別中，如Gypsy－like（6.58%）以及Copia－like（6.84%）LTR型（長(zhǎng)末端重復(fù)）轉(zhuǎn)座元件。大約99.11%的轉(zhuǎn)座元件有超過(guò)10%的分化度，這個(gè)結(jié)果暗示了大多數(shù)的桑樹(shù)轉(zhuǎn)座元件是比較古老的（附圖4）。

1.3 基因預(yù)測(cè)及功能注釋

結(jié)合5個(gè)組織共計(jì)21GB的轉(zhuǎn)錄組數(shù)據(jù)以及5 833個(gè)唯一的表達(dá)序列標(biāo)簽對(duì)基因模型進(jìn)行預(yù)測(cè)及驗(yàn)證，我們?cè)谏?shù)基因組中鑒定到27 085個(gè)高置信度、具有完整基因結(jié)構(gòu)的蛋白編碼基因（附方法及附表6）。在這27 085個(gè)基因中，從頭基因預(yù)測(cè)支持其中99.93%的基因，轉(zhuǎn)錄組測(cè)序及表達(dá)序列標(biāo)簽支持其中58.38%（15 811）的基因，基于同源方法搜索支持其中69.96%（18 943）的基因。這三種方法同時(shí)支持的有超過(guò)一半（52.19%）的基因。包括2 253個(gè)通過(guò)轉(zhuǎn)錄組測(cè)序以及表達(dá)序列標(biāo)簽注釋得到的基因在內(nèi)，我們共預(yù)測(cè)到29 338個(gè)基因（附表7）。這些基因的平均mRNA長(zhǎng)度為2 849bp，平均編碼基因長(zhǎng)度為1 156bp，平均有4.6個(gè)外顯子（附表8）。在這些基因中，轉(zhuǎn)錄組測(cè)序支持其中60.8%的基因，76.92%（22566/29338）的基因在功能數(shù)據(jù)庫(kù)中具有同源的靶基因，如NCBI Nr（NCBI非冗余蛋白質(zhì)數(shù)據(jù)庫(kù)），Swissprot，InterPro，KEGG（京都基因和基因組百科全書(shū)）和COG（直系同源基因蔟）等數(shù)據(jù)庫(kù)（附表9）。

基于轉(zhuǎn)錄組數(shù)據(jù)，我們計(jì)算了組織特異性指數(shù)τ以篩選組織特異性基因和看家基因。我們發(fā)現(xiàn)241、213、285、360和404個(gè)在根、皮、冬芽、雄花和葉中特異表達(dá)基因。同時(shí)，我們發(fā)現(xiàn)1 805個(gè)基因在這5個(gè)組織當(dāng)中持續(xù)表達(dá)，其中包括116個(gè)編碼核糖體蛋白和26個(gè)編碼轉(zhuǎn)錄起始因子的基因（附圖5）。

1.4 基因組進(jìn)化

通過(guò)比較桑樹(shù)與薔薇目（大麻［9］、蘋(píng)果［10］和草莓［11］）基因組序列，可以為了解這一重要種群DNA水平上的差異性提供深刻見(jiàn)解。通過(guò)對(duì)桑樹(shù)和12個(gè)已測(cè)序的植物（圖2）單拷貝基因的系統(tǒng)發(fā)生分析表明?？婆c薔薇科的親緣關(guān)系最近［12－13］。這個(gè)結(jié)果同時(shí)表明桑樹(shù)和大麻兩個(gè)物種的分化是在63.5百萬(wàn)年前，桑樹(shù)和蘋(píng)果兩個(gè)物種的分化是在88.2百萬(wàn)年前，桑樹(shù)和苜蓿的分化則是在101.6百萬(wàn)年前［14］。同義替換分布說(shuō)明桑樹(shù)（?？疲┡c大麻的分化時(shí)間比桑樹(shù)與薔薇科中蘋(píng)果和草莓的分化時(shí)間晚（圖3）。

我們用幾種植物中不同的基因群來(lái)構(gòu)建了3個(gè)系統(tǒng)發(fā)生樹(shù)。首先，桑樹(shù)中的單拷貝基因以及這些單拷貝基因在其他物種中最佳匹配的基因用來(lái)構(gòu)建系統(tǒng)發(fā)生樹(shù)（圖4a）。其次，用genewise預(yù)測(cè)的單拷貝基因構(gòu)建系統(tǒng)發(fā)生樹(shù)（圖4b）。再次，用不同基因組中共線性位置中的最佳匹配的序列構(gòu)建系統(tǒng)發(fā)生樹(shù)。在所有構(gòu)建的系統(tǒng)發(fā)生樹(shù)中，桑樹(shù)分支比其他物種的分支長(zhǎng)，這個(gè)結(jié)果表明桑樹(shù)比其他薔薇目物種的進(jìn)化速度約快3倍。

在沒(méi)有可用的遺傳圖譜的情況下，我們使用基于計(jì)算機(jī)的基因標(biāo)記及基因組組裝方法來(lái)比較桑樹(shù)與草莓基因組的共線性及進(jìn)化關(guān)系［15］。桑樹(shù)和草莓保守共線性區(qū)域基因的分布密度是依據(jù)滑動(dòng)窗口的方法計(jì)算并且以熱度圖的方式進(jìn)行可視化（圖5，附數(shù)據(jù)集1）。

最初通過(guò)對(duì)葡萄基因組的分析發(fā)現(xiàn)了雙子葉植物普遍存在三倍化事件［16］，這里我們將桑樹(shù)scaffolds與桑樹(shù)和葡萄最佳匹配的染色體區(qū)域進(jìn)行比對(duì)可以發(fā)現(xiàn)另外兩個(gè)不太明顯的同源區(qū)域，這個(gè)結(jié)果暗示桑樹(shù)同樣經(jīng)歷了雙子葉植物常見(jiàn)的三倍化事件。草莓和大麻中也發(fā)現(xiàn)一個(gè)主要的和兩個(gè)次要的與葡萄具有同源性的區(qū)域（圖6b，c）。草莓中的這個(gè)結(jié)果與早期報(bào)道草莓中無(wú)古多倍化事件的結(jié)果相反。基因組內(nèi)和基因組間共線性區(qū)域內(nèi)同源基因的同義核苷酸替換率分布進(jìn)一步證明桑樹(shù)、草莓和大麻經(jīng)歷了泛雙子葉植物的六倍化事件（圖6d，c）。

圖2 13種植物的系統(tǒng)發(fā)生關(guān)系分析

圖3 Ks分布圖

圖4 川桑與其他植物的系統(tǒng)發(fā)生樹(shù)

圖5 基于草莓基因組，利用電腦模擬了川桑的基因模型

圖6 物種之間的點(diǎn)圖和Ks分布

1.5 分化選擇

多種薔薇目栽培種中不同的形態(tài)和植物化學(xué)性質(zhì)可能反映了同源基因的分化選擇。通過(guò)對(duì)ω（非同義核苷酸替換率（Ka）對(duì)同義核苷酸替換率（Ks）的比值（Ka/Ks））和Ks值的回歸分析，在桑樹(shù)－大麻、桑樹(shù)－草莓、桑樹(shù)－蘋(píng)果和桑樹(shù)－苜蓿之間鑒定出307、338、353和197個(gè)分化性選擇基因?qū)Γǜ綌?shù)據(jù)集2），其ω顯著高于平均值。有趣的是，用更加嚴(yán)格的Fisher精確性檢驗(yàn)進(jìn)行分析發(fā)現(xiàn)，桑樹(shù)－大麻的222對(duì)分化性選擇基因（附圖6和附表10）主要集中在老齡化和壓力應(yīng)答相關(guān)方面，這可能與植物間生命期長(zhǎng)短不同有關(guān)。桑樹(shù)－草莓，桑樹(shù)－蘋(píng)果的228和258對(duì)分化性選擇的同源基因可能與不同的功能相關(guān)，例如Morus000754（桑樹(shù)）－MDP0000252168（蘋(píng)果）和Morus009486（桑樹(shù)）－MDP0000290357（蘋(píng)果）涉及角皮質(zhì)的合成過(guò)程，可能與蘋(píng)果表皮厚度有關(guān)（盡管表皮生物合成機(jī)制尚不清楚）。在桑樹(shù)－薔薇目（蘋(píng)果，草莓）中受到分化性選擇特別突出的基因?qū)κ桥c質(zhì)體構(gòu)成有關(guān)的（附數(shù)據(jù)集3和4），推測(cè)核酮糖二磷酸羥化酶和許多質(zhì)體基因受到正向的分化選擇。

1.6 抗性基因

桑樹(shù)基因組有142個(gè)具有核苷酸結(jié)合位點(diǎn)（NBS）的抗性基因，占桑樹(shù)所有基因的0.53%，與擬南芥（0.52%）和草莓（0.58%）中所占比例相當(dāng)，但是低于楊樹(shù)（0.86%）和蘋(píng)果（1.49%）（附數(shù)據(jù)集5，附表11）。所有的 R 基因分成6個(gè)類(lèi)群：TIR－NBS－LRR、CC－NBS－LRR、NBSLRR、NBS、CC－NBS以及 TIR－NBS，其中數(shù)目最多的類(lèi)群為CC－NBS－LRR，包括46個(gè)基因。桑樹(shù)基因組中含有127個(gè)半胱氨酸蛋白酶（CP；0.47%）以及129個(gè)天冬氨酸蛋白酶（AP；0.48%）的編碼基因，與蘋(píng)果（0.59%，0.37%）和草莓（0.49%，0.53%）基因組中其所占比例相當(dāng)（附數(shù)據(jù)集6，7和附表12）。其中很突出的一點(diǎn)是有13個(gè)CP和4個(gè)AP基因在桑樹(shù)的乳汁管中表達(dá)（附表13）。有趣的是，四個(gè)AP基因中的一個(gè)（Morus008067）相對(duì)蘋(píng)果中其同源基因（MDP0000201076；附數(shù)據(jù)集2）而言受到分化性選擇。

1.7 蛋白酶抑制子基因

為了減輕昆蟲(chóng)的侵襲損害，植物進(jìn)化出了一種通過(guò)表達(dá)植物蛋白酶抑制子（PIs）干涉昆蟲(chóng)消化系統(tǒng)的防御機(jī)制。在已知的PI序列和它們的保守結(jié)構(gòu)域的基礎(chǔ)上，我們?cè)谏?shù)中鑒定了79個(gè)PIs（附表4）。在桑樹(shù)基因組中注釋出22個(gè)C1家族的半胱氨酸蛋白酶抑制劑基因和19個(gè)A1/C1家族的絲氨酸蛋白酶抑制劑基因，占被鑒定的蛋白酶抑制劑基因的一半。

1.8 家蠶中鑒定到的桑樹(shù)miRNAs

家蠶對(duì)桑葉季節(jié)性生長(zhǎng)的適應(yīng)可能涉及跨界的分子信號(hào)。通過(guò)桑樹(shù)基因組與多個(gè)植物small RNA數(shù)據(jù)庫(kù)的比對(duì)，我們預(yù)測(cè)得到311個(gè)小核RNAs和233個(gè)miRNAs（附表15）。其中有5個(gè)在桑樹(shù)基因組中存在，但在家蠶基因組中不存在的miRNAs，出現(xiàn)在家蠶幼蟲(chóng)血淋巴（2），前中部絲腺（2）和后部絲腺（1）（附表16）。采用另一批家蠶血淋巴對(duì)其進(jìn)行重復(fù)測(cè)序，驗(yàn)證了家蠶血淋巴中存在桑樹(shù)來(lái)源的miRNAs。

2 討 論

早期的研究推測(cè)桑樹(shù)的基礎(chǔ)染色體數(shù)目是1419，盡管之后對(duì)兩種印度桑的細(xì)胞學(xué)實(shí)驗(yàn)表明桑樹(shù)的染色體基數(shù)有可能是7［20］，但是人們一直認(rèn)為桑樹(shù)的染色體基數(shù)是14［19］，并在大量文章引用該觀點(diǎn)。該屬植物的多倍性反映在其染色體數(shù)目的廣泛性上：如川桑（M.notabilis）14條［21］，印度桑（M.indica）、白桑（M.alba）28條，山桑（M.bombycis）42條，黑桑（M.nigra）308條［22］。由于多倍體基因組的復(fù)雜度高，所以選用14條染色體的川桑進(jìn)行全基因組測(cè)序分析。為了證實(shí)川桑染色體數(shù)目，對(duì)處于有絲分裂中期的頂芽細(xì)胞進(jìn)行細(xì)胞學(xué)分析，確認(rèn)了川桑染色體數(shù)目為14條。染色體核型分析表明，川桑的14條染色體能清晰地組成7對(duì)，該結(jié)果支持印度桑關(guān)于桑樹(shù)染色體基數(shù)為7的觀點(diǎn)［20］。

系統(tǒng)進(jìn)化分析法表明桑樹(shù)與其它薔薇目的物種聚為一個(gè)進(jìn)化支。?？瓢凑諔T例認(rèn)為是屬于蕁麻目，而蕁麻目被認(rèn)為是與薔薇目進(jìn)化關(guān)系最近的一個(gè)目。然而最近的一個(gè)研究認(rèn)為榆科、大麻科、?？坪褪n麻科屬于單一的進(jìn)化支［23］，將其命名為urticalean rosids［24］。隨后被子植物種系發(fā)生群III將?？茪w屬為薔薇目［13］，本研究也支持這種重新分類(lèi)的結(jié)果。

桑樹(shù)在核酸水平上進(jìn)化速度很快，其快速進(jìn)化的基因可能使桑樹(shù)不僅對(duì)本地環(huán)境有很好的適應(yīng)性，還能促進(jìn)其傳播到歐洲，非洲和美國(guó)。與桑樹(shù)快速變化的核酸水平相比，薔薇目的倍性進(jìn)化則非常保守。桑樹(shù)，草莓，大麻，番木瓜和葡萄都經(jīng)歷過(guò)最近的泛雙子葉植物六倍化事件。桑樹(shù)的染色體數(shù)目高達(dá)308條（44×）［22］，草莓高達(dá)70條，廣泛存在的高倍數(shù)性揭示了這些譜系能接受新一輪多倍化帶來(lái)的益處。

桑樹(shù)是一種多年生的木本植物，對(duì)桑樹(shù)持續(xù)剪伐不僅能采集桑葉養(yǎng)蠶，還能提高桑葉產(chǎn)量，但是剪伐會(huì)增加桑樹(shù)受到害蟲(chóng)侵?jǐn)_和病原菌侵染的風(fēng)險(xiǎn)，因此，桑樹(shù)需要強(qiáng)大的防御系統(tǒng)抵擋剪伐帶來(lái)的生物脅迫。R基因編碼的蛋白能識(shí)別病原菌的效應(yīng)子，比如其相應(yīng)的無(wú)病毒的基因產(chǎn)物［25］。關(guān)于植物R基因的研究大多數(shù)集中在包含NBS結(jié)構(gòu)域的R基因上［26］。在桑樹(shù)基因組中我們共鑒定到142個(gè)含NBS結(jié)構(gòu)域的R基因。桑樹(shù)是一種含有乳汁的植物，乳汁中的蛋白組分，比如幾丁質(zhì)酶樣蛋白涉及抵御微生物或植食性昆蟲(chóng)的侵?jǐn)_［27－29］。番木瓜乳汁管中的半胱氨酸蛋白酶和豬籠草的天門(mén)冬氨酸蛋白酶對(duì)植食性昆蟲(chóng)均有毒性［30－31］。通過(guò)全基因組分析發(fā)現(xiàn)桑樹(shù)有127個(gè)半胱氨酸蛋白酶編碼基因和129個(gè)天冬氨酸蛋白酶編碼基因，這些基因的功能研究將會(huì)擴(kuò)展我們對(duì)桑樹(shù)防御機(jī)制的理解。

植物防御機(jī)制能干擾昆蟲(chóng)的消化系統(tǒng)，關(guān)于寡食性家蠶是如何繞過(guò)桑樹(shù)的防御機(jī)制，這個(gè)問(wèn)題至今為止仍不清楚。植物蛋白酶抑制劑能降低植食性昆蟲(chóng)中腸消化酶的活性，從而導(dǎo)致昆蟲(chóng)發(fā)育嚴(yán)重畸形，死亡，生殖率降低［32－33］。早期研究報(bào)道了植物產(chǎn)生的多結(jié)構(gòu)域、多聚體結(jié)構(gòu)的蛋白酶抑制劑對(duì)斜紋夜蛾具有抗?fàn)I養(yǎng)作用［34］。昆蟲(chóng)可以通過(guò)誘導(dǎo)產(chǎn)生一些對(duì)蛋白酶抑制劑不敏感的蛋白酶以及用體內(nèi)特殊的蛋白酶降解植物的蛋白酶抑制劑等方法來(lái)避開(kāi)植物蛋白酶抑制劑的作用［35－36］。一種對(duì)十字花科有害的著名鱗翅目昆蟲(chóng)—小菜蛾能使芥末的胰蛋白酶抑制劑2失活，從而打破宿主植物的抵御［37］。受益于家蠶和桑樹(shù)的基因組序列，對(duì)其轉(zhuǎn)錄組的比較分析能夠推進(jìn)我們對(duì)植物－植食性昆蟲(chóng)相互適應(yīng)關(guān)系的理解。

在家蠶三個(gè)組織的小RNA測(cè)序數(shù)據(jù)中發(fā)現(xiàn)五個(gè)來(lái)源于桑樹(shù)而非家蠶的miRNA，其中的MIR156在水稻營(yíng)養(yǎng)生長(zhǎng)期的老葉中高量表達(dá)，該miRNA在植物從幼態(tài)向成熟態(tài)的轉(zhuǎn)變中起主要作用［38－40］。同時(shí)我們注意到水稻的MIR168a能轉(zhuǎn)運(yùn)到人體內(nèi)，并調(diào)節(jié)低密度脂蛋白受體銜接蛋白1的表達(dá)水平［41］，那么存在家蠶絲腺中的桑樹(shù)MIR156是否向家蠶發(fā)出桑葉衰老的信號(hào)，促進(jìn)家蠶吐絲成繭，亦或是家蠶組織特異存在的其它桑樹(shù)miRNA在家蠶的發(fā)育中是否起作用，這些問(wèn)題至今仍不清楚。

綜上所訴，基因組信息對(duì)現(xiàn)代桑樹(shù)遺傳研究是十分重要的資源。桑樹(shù)基因組的特點(diǎn)，如基因家族，片段重復(fù)，共線性區(qū)域等不僅豐富了植物比較基因組可用的數(shù)據(jù)，并且促進(jìn)了與?？泼芮邢嚓P(guān)物種靶基因的鑒定?；诨蚪M序列發(fā)展的遺傳標(biāo)記可以用于構(gòu)建遺傳圖譜，定位克隆，品系鑒定以及標(biāo)記輔助的篩選。這些分子工具和基因組技術(shù)的應(yīng)用將會(huì)有效促進(jìn)農(nóng)業(yè)的發(fā)展。桑樹(shù)和家蠶作為植物－植食性動(dòng)物相互作用的模式系統(tǒng)，二者的基因組序列對(duì)于更深層次解析普遍存在大多數(shù)陸生環(huán)境中的生物學(xué)合伙關(guān)系提供了一個(gè)獨(dú)一無(wú)二的機(jī)會(huì)。

3 研究方法

3.1 川桑核型分析

川桑嫩葉用2mM 8－羥基喹啉在室溫下處理3h，之后用3：1的甲醇：冰醋酸在4℃固定2h。固定好的葉片用1/15MKCl溶液處理30min，隨后用2.5%（W/V）的纖維素酶和2.5%（W/V）的果膠酶（YaKult Co.，日本）在37℃消化1.5h。消化后的葉片用ddH2O處理10min，之后再用3：1的甲醇：冰醋酸在室溫下固定30min。再固定的葉片打碎后滴兩滴懸浮液在載玻片上，室溫下用吉姆薩染液染色6小時(shí)，隨后在顯微鏡下進(jìn)行觀察（Olympus，日本）。

3.2 DNA和RNA制備

測(cè)序用桑樹(shù)品種為野生桑種川桑，其含有14條染色體。川桑冬芽的基因組DNA用CTAB法進(jìn)行提取并用于測(cè)序文庫(kù)的構(gòu)建。五個(gè)組織（根、一年生皮、冬芽、雄花和葉）的總RNA參照Wan和Wilhins等的方法進(jìn)行提?。?2］，并用不含RNA酶的DNA酶Ⅰ（New England BioLabs）在37℃下消化30min以去除剩余的DNA。用含有oligo（dT）的珠子分離含有poly（A）的mRNA。用隨機(jī)的六聚引物和反轉(zhuǎn)錄酶（Invitrogen）合成cDNA第一鏈。用DNA聚合酶（New England BioLabs）和RNA酶H（Invitrogen）合成cDNA第二鏈。

3.3 基因組測(cè)序

桑樹(shù)基因組測(cè)序使用全基因組鳥(niǎo)槍法。測(cè)序文庫(kù)根據(jù)制造商（Illumina，圣地亞哥，加拿大）的說(shuō)明進(jìn)行制備。對(duì)于短插入片段的DNA文庫(kù)，將5μg基因組DNA用壓縮氮?dú)鈬婌F法進(jìn)行片段化處理。片段化的DNA末端用堿基A進(jìn)行鈍化處理，然后將3’末端含有一個(gè)堿基T的DNA接頭（Illumina）與DNA片段進(jìn)行連接。之后用2%瓊脂糖凝膠將連接產(chǎn)物進(jìn)行分離，通過(guò)切膠純化各個(gè)插入大小的DNA片段。對(duì)于長(zhǎng)的mate－paired文庫(kù)（大于等于2Kb），將10～30μg基因組DNA用壓縮氮?dú)鈬婌F法進(jìn)行片段化處理，然后用生物素標(biāo)記的dNTPs對(duì)片段進(jìn)行處理，通過(guò)凝膠選擇2Kb，5Kb和10Kb的主要條帶。之后DNA片段通過(guò)自連而環(huán)化。將DNA片段的兩端進(jìn)行結(jié)合并用核酸外切酶對(duì)線性DNA片段進(jìn)行消化。消化之后環(huán)化的DNA片段重新進(jìn)行片段化，隨后利用生物素與抗生蛋白鏈菌素的相互作用，用磁珠富集含有之前結(jié)合的末端的DNA，之后將這些末端進(jìn)行鈍化并加上含有堿基A的接頭。我們根據(jù)制造商的說(shuō)明進(jìn)行了如下的paired－end（PE）測(cè)序流程：成蔟、模板雜交、等溫?cái)U(kuò)增、線性化、封閉以及測(cè)序引物的變性和雜交。測(cè)序完成后用base－calling程序（SolexaPipeling－0.3 ）從最初的熒光圖像中獲得序列。

3.4 基因組組裝

在進(jìn)行從頭組裝之前，我們通過(guò)以下五個(gè)步驟對(duì)低質(zhì)量的數(shù)據(jù)進(jìn)行過(guò)濾：（1）根據(jù)Hiseq2000的測(cè)序質(zhì)量報(bào)告去除5’和3’末端的低質(zhì)量的堿基；（2）去除含有N大于10%的序列；（3）去除低質(zhì)量堿基（Q＜8）大于50%的序列；（4）去除被接頭污染的序列；（5）去除在文庫(kù)構(gòu)建過(guò)程中因PCR引起的重復(fù)的序列。SOAPdenovo是由深圳華大基因研發(fā)的一種基因組組裝軟件，該軟件采用基于圖論的算法和分步組裝的策略［43］。我們首先用49－kmers對(duì)小插入片段的文庫(kù)（＜1kb）進(jìn)行組裝。然后我們用41－kmers將所有過(guò)濾得到的序列重新比對(duì)到組裝好的重疊群上并且將比對(duì)上的序列編譯為可用的重疊群。根據(jù)paired－end信息，我們共采用了七步，從170bp的插入文庫(kù)到20kb的插入文庫(kù)，將重疊群組裝成scaffold。為了填補(bǔ)scaffold中的空隙，我們收集了paired－end的序列，如果其中一端比對(duì)到一個(gè)重疊群上而另一端位于空隙中，則重復(fù)進(jìn)行一次局部的組裝。

3.5 轉(zhuǎn)座元件（TE）和重復(fù)序列

為了預(yù)測(cè)桑樹(shù)基因組中的TE，我們首先用RepeatModeler（version 1.0.3，http：//www.repeatmasker.org/RepeatModeler.html），RepeatScout［44］（version 1.0.5，http：//bix.ucsd.edu/repeatscout/）和Piler［45］（version 1.0，http：//www.drive5.com/piler/）構(gòu)建了一個(gè) TE庫(kù)，之后用 RepeatMasker［46］（version 3.2.9，http：//www.repeatmasker.org/）對(duì)桑樹(shù)基因組中的TE進(jìn)行了一個(gè)從頭預(yù)測(cè)。我們同時(shí)也使用RepeatMasker和ProteinMask（version 3.0）對(duì)已知類(lèi)型的TE進(jìn)行了預(yù)測(cè)，在預(yù)測(cè)過(guò)程中我們使用的TE庫(kù)包括Repbase［47］（version 15.02，http：//www.girinst.org/repbase/）和TIGR中的雙子葉植物 TE庫(kù)［48］（version 3.0，ht－tp：//plantta.jcvi.org）。串聯(lián)重復(fù)序列使用Tandem Repeats Finder（TRF，version 4.04，http：//tandem.bu.edu/trf/trf.html）來(lái)進(jìn)行預(yù)測(cè)。簡(jiǎn)單重復(fù)序列、微衛(wèi)星序列和低復(fù)雜度重復(fù)序列使用RepeatMasker進(jìn)行預(yù)測(cè)（使用“－noint”選項(xiàng)）［49］。將桑樹(shù)基因組中已進(jìn)行分類(lèi)的TE家族與Repbase（v15.02）中的共有序列進(jìn)行比對(duì)以確定這些TE的序列分化度。

3.6 基因預(yù)測(cè)和注釋

桑樹(shù)基因預(yù)測(cè)共使用了三種方法：基于同源性比較的方法、從頭預(yù)測(cè)方法和基于表達(dá)序列標(biāo)簽（EST）/轉(zhuǎn)錄組數(shù)據(jù)的方法。其中高置信度的基因是用前兩種方法預(yù)測(cè)得到。為了注釋桑樹(shù)中的蛋白編碼基因，我們將其27085個(gè)高置信度基因的核苷酸序列與NCBI、KEGG、COG和Swissport等數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，比對(duì)的閾值設(shè)為1e－5。蛋白質(zhì)結(jié)構(gòu)域和功能預(yù)測(cè)在Iprscan（v4.4.1）上進(jìn)行注釋。

3.7 RNA和EST測(cè)序

cDNA文庫(kù)的制備和測(cè)序參照Illumina的操作流程。用TopHat（v1.3.3）將RNA測(cè)序得到的序列與桑樹(shù)基因組進(jìn)行比對(duì)。通過(guò)計(jì)算每百萬(wàn)條reads中覆蓋到每kb上的reads數(shù)來(lái)衡量五個(gè)組織中的基因表達(dá)水平，并且通過(guò)計(jì)算組織特異性的指數(shù)τ來(lái)鑒定組織特異性表達(dá)的基因。對(duì)于EST測(cè)序，將相同的五個(gè)組織的RNA混合并用Creator SMART cDNA Kit（Clontech）進(jìn)行cDNA的合成。用Trimmer－Director kit（Evrogen）構(gòu)建標(biāo)準(zhǔn)化的cDNA文庫(kù)。從該文庫(kù)中隨機(jī)挑選10 000個(gè)克隆用ABI3730（Applied Biosystem）進(jìn)行測(cè)序。

3.8 非編碼RNA基因

桑樹(shù)基因組中的tRNA使用tRNAscan－SE（v1.23）軟件進(jìn)行預(yù)測(cè)（使用“eukaryotes”選項(xiàng)）［50］。用BLASTN將桑樹(shù)基因組與植物rRNA序列進(jìn)行比對(duì)（閾值為1e－5），序列一致性大于85%并且高分值片段對(duì)的長(zhǎng)度大于50bp的結(jié)果認(rèn)為是桑樹(shù)的rRNA。對(duì)于microRNA和snRNA，首先用BLASTN將桑樹(shù)基因組與Rfam數(shù)據(jù)庫(kù)（v9.1）進(jìn)行比對(duì)（閾值為1），比對(duì)結(jié)果進(jìn)一步用INFERNAL軟件進(jìn)行分析，該軟件根據(jù)RNA的結(jié)構(gòu)和序列相似性對(duì)microRNA和snRNA進(jìn)行預(yù)測(cè)。

3.9 模擬基因染色

我們通過(guò)BLASTP比對(duì)（閾值為1e－5）鑒定桑樹(shù)和草莓之間的直系同源基因?qū)?，得到的結(jié)果進(jìn)一步用Mcscan分析以鑒定兩個(gè)物種間的共線性區(qū)域。通過(guò)使用Genome Zipper軟件，根據(jù)草莓的連鎖群，將桑樹(shù)中與之對(duì)用的共線性區(qū)域的scaffold進(jìn)行連接。計(jì)算每500kb滑動(dòng)窗口內(nèi)的基因密度和對(duì)應(yīng)的直系同源基因的密度并將結(jié)果分布做熱度圖。

3.10 在家蠶的組織中鑒定桑樹(shù)的miRNAs

用mirVana PARIS kit（Ambion，美國(guó)）從12ml的家蠶血淋巴（取自五齡5天幼蟲(chóng)）中提取小RNA。提取的小RNA測(cè)序參照Liu等人的操作流程進(jìn)行［51］。前中部和后部絲腺的小RNA數(shù)據(jù)從http：//www.ncbi.nlm.nih.gov/gds？term＝GSE17965上下載。這三個(gè)家蠶組織的小RNA序列通過(guò)BLASTN與桑樹(shù)預(yù)測(cè)的miRNAs進(jìn)行比對(duì)（完全匹配）以鑒定家蠶組織中存在的桑樹(shù)來(lái)源的miRNAs。

3.11 系統(tǒng)發(fā)生樹(shù)的構(gòu)建和物種形成時(shí)間的確定

基于最大似然法，我們利用13個(gè)植物物種的單拷貝基因構(gòu)建了系統(tǒng)發(fā)生樹(shù)。通過(guò)BLAST比對(duì)（閾值為1e－10）確定物種間的直系同源基因?qū)Γ鼈冎g的同義替換率（Ks）［52］使用PAML軟件包中的yn00程序進(jìn)行計(jì)算［53］?；贙s值的物種形成時(shí)間根據(jù)以下公式計(jì)算：T＝ Ks/2λ（其中λ＝ 6.1×10－9）［54］。根據(jù)非同義替換率（Ka）和 Ks之間的回歸分析（95%預(yù)測(cè)區(qū)間范圍）［55］推測(cè)桑樹(shù)和其他四種植物間可能受到正選擇的直系同源基因?qū)?。ω值大于預(yù)測(cè)區(qū)間上限的基因?qū)φJ(rèn)為是受到正選擇的。用BLAST2GO［56］對(duì)具有高ω值的基因?qū)Φ幕虮倔w論類(lèi)群進(jìn)行預(yù)測(cè)，其結(jié)果用Fisher精確性檢驗(yàn)進(jìn)行驗(yàn)證（P值＜0.05）。

3.12 推測(cè)基因共線性

我們用MCSCAN57（一種多染色體比對(duì)工具）和COLINEARSCAN58（一種成對(duì)染色體比對(duì)工具）相結(jié)合來(lái)推測(cè)共線性基因，這些共線性基因用于系統(tǒng)發(fā)生和進(jìn)化分析。

3.13 推測(cè)進(jìn)化事件的時(shí)間

我們利用物種間的共線性基因和物種內(nèi)高度可信的同源基因來(lái)推測(cè)進(jìn)化事件。例如，桑樹(shù)scaffold之間的共線性基因可能來(lái)源于祖先的多倍化事件（如果存在的話），而桑樹(shù)和葡萄之間的共線性基因可能來(lái)源于兩個(gè)物種的分化。Ks使用PAML［53］中的Nei－Gojobori方法［52］計(jì)算。通過(guò)描繪Ks值的分布來(lái)推測(cè)進(jìn)化事件的相對(duì)時(shí)間。

3.14 同源基因點(diǎn)圖

在這一部分分析中我們使用之前介紹的預(yù)測(cè)基因集和Genewise［59］預(yù)測(cè)的一個(gè)基因集。葡萄、蘋(píng)果、草莓和大麻的基因組序列和注釋下載自在線數(shù)據(jù)庫(kù)，所有的數(shù)據(jù)使用2012年10月之前的最新版本。在比較有可用的擬染色體的基因組時(shí)，我們用BLASTP進(jìn)行蛋白－蛋白之間的比對(duì)，從而找到可能的同源基因，得到的結(jié)果用于作點(diǎn)圖；作圖過(guò)程中以基因在染色體上的順序?yàn)樽鴺?biāo)。當(dāng)比較沒(méi)有可用的擬染色體的基因組時(shí)（例如桑樹(shù)和大麻），我們選擇一個(gè)有可用的擬染色體物種（例如葡萄）的編碼基因序列，用BLASTN與桑樹(shù)和大麻的基因組進(jìn)行比對(duì)，比對(duì)的結(jié)果用于作點(diǎn)圖。為了推測(cè)基因組重復(fù)事件，沒(méi)有錨定的scaffold根據(jù)它們?cè)谄咸鸦蚪M上的最好匹配區(qū)域連接成假定的擬染色體，桑樹(shù)和大麻假定的擬染色體就是用這種方法構(gòu)建。在點(diǎn)圖中，相應(yīng)的葡萄的基因區(qū)域?qū)?duì)應(yīng)兩個(gè)匹配的區(qū)域。

3.15 本研究中所使用的數(shù)據(jù)

本研究中所使用的基因組數(shù)據(jù)從以下網(wǎng)址中下載，并列出了其在NCBI上的登錄號(hào)。擬南 芥 （TAIR9）：ftp：//ftp.arabidopsis.org/Genes/TAIR9＿genome＿release/，GCA ＿000001735.1。大麻：http：//genome.ccbr.utoronto.ca/downloads.html，GCA＿000230575.1。番木 瓜 （version 1th）：ftp：//ftp.jgi－psf.org/pub/compgen/phytozome/v5.0/Cpapaya/，GCA＿000150535.1。黃瓜（version 1th）：http：//cucumber.genomics.org.cn/page/cucumber/download.jsp，GCA＿000004075.1。草莓（version 1.1）：http：//www.rosaceae.org/species/fragaria/fragaria＿vesca/genome＿v1.1，GCA＿000184155.1。大豆（version 1.0）：ftp：//ftp.jgi－psf.org/pub/compgen/phytozome/v5.0/Gmax/，GCA＿000004515.1。蘋(píng)果（version 1.0）：http：//genomics.research.iasma.it/index.html，GCA＿000148765.2。苜蓿：ftp：//ftp.jgi－psf.org/pub/compgen/phytozome/v8.0/Mtruncatula/，GCA＿000219495.1。楊樹(shù) （version 5.0）：ftp：//ftp.jgi－psf.org/pub/compgen/phytozome/v5.0/Ptrichocarpa/，GCA ＿000002775.1。桃：ftp：//ftp.jgi－psf.org/pub/compgen/phytozome/v8.0/Ppersica/，GCA＿000346465.1。白梨：http：//peargenome.njau.edu.cn：8004/default.asp？d＝1＆m＝1，GCA＿000315295.1?？煽桑╲ersion 1.0）：http：//cocoagendb.cirad.fr/gbrowse/download.html，GCA＿000403535.1。葡萄：http：//www.genoscope.cns.fr/externe/Download/Projets/Projet＿M(jìn)L/data/12X/，GCA＿000003745.2。

論 文 鏈 接：http：//www.nature.com/ncomms/2013/130919/ncomms3445/full/ncomms3445.html.