桑樹（川桑Morus notabilis）全基因組測序＊

桑樹基因組研究團隊

（家蠶基因組生物學國家重點實驗室，重慶 400716）

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

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

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

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

1 結 果

1.1 基因組測序和組裝

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

圖1 川桑（M.notabilis）染色體細胞學分析

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

1.2 重復序列

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

1.3 基因預測及功能注釋

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

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

1.4 基因組進化

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

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

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

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

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

圖3 Ks分布圖

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

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

圖6 物種之間的點圖和Ks分布

1.5 分化選擇

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

1.6 抗性基因

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

1.7 蛋白酶抑制子基因

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

1.8 家蠶中鑒定到的桑樹miRNAs

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

2 討 論

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

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

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

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

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

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

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

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小時，隨后在顯微鏡下進行觀察（Olympus，日本）。

3.2 DNA和RNA制備

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

3.3 基因組測序

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

3.4 基因組組裝

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

3.5 轉座元件（TE）和重復序列

為了預測桑樹基因組中的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/）構建了一個 TE庫，之后用 RepeatMasker［46］（version 3.2.9，http：//www.repeatmasker.org/）對桑樹基因組中的TE進行了一個從頭預測。我們同時也使用RepeatMasker和ProteinMask（version 3.0）對已知類型的TE進行了預測，在預測過程中我們使用的TE庫包括Repbase［47］（version 15.02，http：//www.girinst.org/repbase/）和TIGR中的雙子葉植物 TE庫［48］（version 3.0，ht－tp：//plantta.jcvi.org）。串聯(lián)重復序列使用Tandem Repeats Finder（TRF，version 4.04，http：//tandem.bu.edu/trf/trf.html）來進行預測。簡單重復序列、微衛(wèi)星序列和低復雜度重復序列使用RepeatMasker進行預測（使用“－noint”選項）［49］。將桑樹基因組中已進行分類的TE家族與Repbase（v15.02）中的共有序列進行比對以確定這些TE的序列分化度。

3.6 基因預測和注釋

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

3.7 RNA和EST測序

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

3.8 非編碼RNA基因

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

3.9 模擬基因染色

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

3.10 在家蠶的組織中鑒定桑樹的miRNAs

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

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

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

3.12 推測基因共線性

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

3.13 推測進化事件的時間

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

3.14 同源基因點圖

在這一部分分析中我們使用之前介紹的預測基因集和Genewise［59］預測的一個基因集。葡萄、蘋果、草莓和大麻的基因組序列和注釋下載自在線數(shù)據(jù)庫，所有的數(shù)據(jù)使用2012年10月之前的最新版本。在比較有可用的擬染色體的基因組時，我們用BLASTP進行蛋白－蛋白之間的比對，從而找到可能的同源基因，得到的結果用于作點圖；作圖過程中以基因在染色體上的順序為坐標。當比較沒有可用的擬染色體的基因組時（例如桑樹和大麻），我們選擇一個有可用的擬染色體物種（例如葡萄）的編碼基因序列，用BLASTN與桑樹和大麻的基因組進行比對，比對的結果用于作點圖。為了推測基因組重復事件，沒有錨定的scaffold根據(jù)它們在葡萄基因組上的最好匹配區(qū)域連接成假定的擬染色體，桑樹和大麻假定的擬染色體就是用這種方法構建。在點圖中，相應的葡萄的基因區(qū)域?qū)獌蓚€匹配的區(qū)域。

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

本研究中所使用的基因組數(shù)據(jù)從以下網(wǎng)址中下載，并列出了其在NCBI上的登錄號。擬南 芥 （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。蘋果（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。楊樹 （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。可可（version 1.0）：http：//cocoagendb.cirad.fr/gbrowse/download.html，GCA＿000403535.1。葡萄：http：//www.genoscope.cns.fr/externe/Download/Projets/Projet＿ML/data/12X/，GCA＿000003745.2。

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