小地老虎不同組織表達差異揭示其翅發育機制

鄭偉剛 靳明輝 Swapan Chakrabarty 肖彬 蕭玉濤 袁海濱

（1. 吉林農業大學農學院，長春 130118；2. 中國農業科學院農業基因組研究所，深圳 518120）

小地老虎（Agrotis ipsilon）、棉鈴蟲（Helicoverpa armigera）、帝王蝶、黏蟲（Mythimna separata）、稻飛虱（Nilaparvata lugens）、玉米螟（Ostriniafurnacalis、Ostrinia nubilalis）及東亞飛蝗（Locusta migratoria manilensis（Meyen））等都屬于遷飛性昆蟲［1-6］。遷飛是動物種群對棲息地資源發生季節性或偶然性變化時進化出的一種適應性對策，昆蟲的遷飛行為嚴重影響人類的生產和生活，如危害農作物，造成農作物大量減產、傳播人和牲畜的疾病等。影響昆蟲遷飛的因素眾多，包括生存環境、激素調控、遺傳因素等，而翅發育是昆蟲遷飛的一個重要影響因素。昆蟲翅膀的進化對昆蟲覓食、求偶、避敵、抵御外界不良環境等起到重要作用［7-8］。隨著果蠅等模式昆蟲的深入研究，一系列與翅發育有關的信號通路逐漸被揭示出來，如Wnt、Hedgehog、Decapentaplegic（Dpp）、Hippo和 c-Jun N-terminal Kinase（JNK）等信號通路，其中Wnt通路是昆蟲翅發育的關鍵信號通路。

Wnt信號通路包含3種Wnt信號傳導通路：經典Wnt通路，非經典Wnt /鈣通路以及非經典planar cell polarity通路［9-11］，這 3種途徑都是通過Wnt蛋白配體與Frizzled（Fz）家族受體結合而激活的。經典Wnt通路引起基因轉錄的調節，并且被認為部分地被SPATS1基因負調節［11］。非經典planar cell polarity通路調控細胞骨架，骨架決定細胞的形狀。非經典Wnt /鈣通路，調控細胞內的鈣。Decapentaplegic（Dpp）信號通路中，Dpp是參與果蠅發育的關鍵形態發生素，也是第一個經過驗證的分泌形態發生素［12-13］。早期果蠅胚胎和15個成蟲器的正確形成和發育需要Dpp調控，這些組織將生長發育為成年果蠅中的四肢和其他器官和結構。DPP在調節組織的生長和器官組織大小方面起著重要的作用［12，14］。Hedgehog（Hh）信號通路在控制細胞增殖、組織模式形成、干細胞維持和發育方面發揮著多種作用［15-18］。Hedgehog信號通路是一種向胚胎細胞傳遞信息的信號通路，這些信息是胚胎細胞分化所必需的。胚胎的不同部位具有不同濃度的Hedgehog基因信號蛋白。該通路在成蟲個體中也有作用。Hippo信號通路也被稱為Salvador或者Warts通路，通過調控細胞增殖和凋亡來控制動物的器官大小。該通路的名字來源于它的一個關鍵信號元件——蛋白激酶Hippo（Hpo）［19-20］。該基因的突變導致組織過度生長，或類似“河馬”表型。

小地老虎（Agrotis ipsilon）又名切根蟲、土蠶等，屬于一種世界性的鱗翅目夜蛾科農業害蟲［5，21］，也是典型的遷飛性農業害蟲。由于翅是在成蟲階段發育形成，且屬于昆蟲胸部的組成部分，本試驗對小地老虎幼蟲，成蟲頭部、胸部、腹部和進行轉錄組測序，探究可能決定小地老虎翅發育的關鍵通路及相關基因［22-23］。

1 材料與方法

1.1 材料

小地老虎幼蟲主要來源于中國農業科學院農業基因組研究所農業生態基因組學中心實驗室，小地老虎成蟲主要來源于山東長島野外采集。

1.2 方法

1.2.1 小地老虎總RNA的提取、文庫構建和測序 小地老虎幼蟲6頭作為一個重復，共取18頭做3組幼蟲整蟲轉錄組Ya、Yb、Yc的組織材料，為避免腸道微生物污染，幼蟲整蟲解剖去除圍食膜，再用液氮研磨。小地老虎成蟲每6頭作為一個重復，分別取頭、胸、腹，做頭（Ta、Tb、Tc）、胸（Xa、Xb、Xc）、腹（Fa、Fb、Fc）轉錄組材料。

取烘干研缽，加入少量乙醇，灼燒防止研缽內殘留其他物質。待研缽烘干后，加入液氮，將組織樣品加入研缽中研磨至粉狀，取適量放入裝有1 mL TRIzol的RNAase free 1.5 mL離心管中。通常情況下，1 mL TRIzol Reagent 能夠裂解50-100 mg組織樣品。

向樣品中加入0.2 mL氯仿（每1 mL TRIzol加0.2 mL氯仿）后，用手輕輕搖動15 s左右，使其混合并于常溫孵育2-3 min，然后將樣品于12 000 r/min、4℃條件下進行離心15 min。將離心管傾斜為45°并用移液器將水相吸出，然后移入新的離心管，在此過程中應盡量避免吸入混合相或者有機層。取0.5mL 異丙醇（每1 mL TRIzol對應0.5 mL異丙醇）加入之前移出的水相樣品中，常溫孵育10 min后在4℃、12 000 r/min條件下，離心10 min，倒掉上清液并留下RNA沉淀。加入1 mL 75%酒精（每1 mL的 TRIzol至少加1 mL的75%乙醇），輕微震蕩混勻后在4℃、7 500 r/min條件下離心5 min，倒掉上清液，只留底部沉淀。風干RNA沉淀5-10 min，禁止真空離心RNA沉淀。加入0.02-0.05 mL無核酸酶水，使用移液器反復吸打，然后在55-60℃條件下水浴10-15 min，使RNA沉淀溶解，將得到的RNA溶液于-80℃條件下長期保存。

利用 Nanodrop2000 對根據上述步驟提取出的total RNA，檢測其濃度和純度，利用瓊脂糖凝膠電泳檢測 RNA 完整性，使用Agilent 2100 測定total RNA RIN 值。單次建庫要求RNA 總量 5 μg，濃度≥ 200 μg/mL，OD260/OD2280介于 1.8-2.2 之間。利用Illumina公司的TruSeqTMRNA sample preparation kit完成cDNA文庫的構建，并通過Illumina Hiseq 測序平臺的150 paired-end模式進行轉錄組測序。

1.2.2 轉錄組組裝和ORF預測 對于無參考基因組的轉錄組測序和相關研究，獲得RNA-seq高質量測序數據后，需要將測序數據進行組裝，才能進行后續注釋及其他生物學功能分析。Trinity（http：//trinityrnaseq.sourceforge.net/，版本號 ：v2.6.6）是目前從頭組裝Illumina 短片段序列比較常用的軟件，本試驗在對原始數據進行過濾接頭和低質量數據后，使用Trinity進行從頭組裝。使用TransDecoder（https：//github.com/TransDecoder/TransDecoder，版本號：5.5.0）軟件識別轉錄本序列中候選的編碼區域，結合 cd-hit（http：//weizhongli-lab.org/cd-hit/，版本號：v4.8.1）進行去冗余，得到獨立基因（Unigene）序列。

1.2.3 轉錄組功能注釋 使用Trinity軟件組裝得到轉錄本序列后，要對生物學相關問題進行探究，就要對這些轉錄本進行功能注釋。將拼接所得所有（unigene）序列，使用Interperscan（http：//www.ebi.ac.uk/interpro/，版本號：v5）進行結構域預測，使用 Blastp（https：//www.ncbi.nlm.nih.gov/）分別與KEGG、eggNOG、NR、Trembl、Swissprot、數據庫進行比對獲得相應的功能注釋信息，然后進行KEGG通路富集分析。KEGG（Kyoto Encyclopedia of Genes and Genomes，京都基因和基因組百科全書，http：//www.genome.jp/kegg/）是基因組注釋方面的公共數據庫。KEGG數據庫中豐富的通路信息將有助于我們從系統水平去了解基因的生物學功能，例如代謝途徑、遺傳信息傳遞以及細胞過程等一些復雜的生物功能，這大大提高了該數據庫在實際生產和應用中的價值。根據與KEGG數據庫比對結果，獲得轉錄本對應的KO編號，根據KO編號可以獲得某轉錄本可能參與的具體生物學通路。eggNOG（evolutionary genealogy of genes：Non-supervised Orthologous Groups）數據庫包含了豐富的注釋信息，包含了KEGG等數據庫信息。Swiss-Prot是一個手工注釋的非冗余蛋白序列數據庫。它結合了從科學文獻中提取的信息和生物尿酸評估計算分析，TrEMBL（翻譯的EMBL核苷酸序列數據庫），為那些不在Swiss-Prot中的蛋白質提供自動注釋。NR（NCBI 非冗余蛋白庫）為綜合數據庫，包括了SwissProt、PIR（Protein information resource）、PRF（Protein research foundation）、PDB（Protein data bank）蛋白質數據庫中非冗余的數據以及從GenBank和RefSeq的CDS數據庫中翻譯所得的蛋白質數據。

1.2.4 差異表達基因分析 使用組裝軟件將轉錄組組裝好后，利用RSEM（http：//deweylab.biostat.wisc.edu/rsem/）進行表達量分析。RSEM利用bowtie的比對結果進行表達量統計，首先得到每個樣品比對到每個基因上的read count數目，然后對其進行TPM數值轉換，TPM數值可以反映基因的表達水平。使用R包（edgeR和DESeq）分別對小地老虎轉錄組進行基因表達差異分析，取兩者交集，得到最終差異基因。

2 結果

2.1 原始測序數據質控

測序小地老虎幼蟲，成蟲頭、胸和腹部，各3個樣本，共計12個樣本。小地老虎幼蟲共測數據19.59 G，小地老虎成蟲頭部20.73 G，小地老虎成蟲胸部21.36 G，小地老虎成蟲腹部21.32 G，共計數據量82.99 G，質控去掉Adapter和低于Q20的數據后最終得到78.67 G（表1）。

2.2 轉錄組組裝和ORF預測

使用 Trinity（http：//trinityrnaseq.sourceforge.net/）軟件，對小地老虎轉錄組數據進行組裝。最終組裝得到轉錄組243.4 M，最長序列27.4 K，平均長度0.6 K，N50 0.7 K，denove轉錄組組裝結果GC含量為39.66%（表2）。使用TransDecoder對轉錄組進行ORF 預測，并去冗余后，得到37 744個unigenes。轉錄組組裝長度頻數分布（圖1），圖1-A為轉錄本（Transcripts）組裝長度頻數分布，在200-500范圍內的轉錄本比例最多，為67.59%，其次為在500-1 000范圍內的轉錄本比例為22.03%。圖1-B為Unigenes組裝長度頻數分布，在500-1 000范圍內的轉錄本比例最多，為33.83%，其次為在200-500范圍內的比例為30.29%。

表1 小地老虎轉錄組測序原始數據及質控數據統計

表2 轉錄組denovo組裝結果統計

2.3 小地老虎轉錄組功能注釋

使用Interperscan進行結構域預測，共有20 522個unigenes預測出結構域信息。將拼接所得獨立基因，使用Blastp分別與NR、Sport、Trembl、KEGG數據庫進行比對，比對后取besthit獲得相應的注釋信息。

KEGG 數據庫注釋功能基因23 334個，注釋通路373條，KO 6 161個；eggNOG 數據路注釋功能基因31 556個，Orthologous Groups（OG，同源群）8 211個；Nr 數據庫共注釋功能基因29 971個；Swissprot 數據庫注釋功能基因本21 538個；Trembl數據庫注釋功能基因32 354個。各個數據庫功能注釋Venn圖（圖2），其中5個數據庫能夠共同注釋到的基因有18 345個。5個數據庫分別注釋出了本數據庫特有基因，其中Trembl數據庫注釋出的特有基因最多，為426個，KEGG數據庫注釋出的特有基因最少，有15個（圖2）。

圖1 轉錄組組裝長度頻數分布

2.4 小地老虎轉錄組表達差異分析

分別使用EdgeR.R和DEseq2.R 進行差異分析，并求出差異基因。再將兩種方法得到的差異基因求交集，得到最終的差異基因（圖3）。最終得到小地老虎成蟲頭部樣本與小地老虎胸部樣本轉錄組差異基因4 069個，小地老虎成蟲胸部樣本與小地老虎腹部樣本轉錄組差異基因6 905個，小地老虎成蟲腹部樣本與小地老虎頭部樣本轉錄組差異基因7 777個，小地老虎幼蟲與小地老虎成蟲頭部樣轉錄組差異基因5 768個，小地老虎幼蟲與小地老虎成蟲胸部樣轉錄組差異基因4 836個，小地老虎幼蟲與小地老虎成蟲腹部樣轉錄組差異基因8 917個。小地老虎成蟲頭部與腹部差異基因最多，小地老虎成蟲頭部與胸部差異基因最少。小地老虎幼蟲與成蟲頭、胸、腹三部分相比，小地老虎幼蟲與成蟲腹部的差異基因最多，與胸部的差異基因最少（圖3）。

圖2 轉錄組不同數據庫注釋結果Venn圖

圖3 各組樣本DESeq2與EdegR基因差異表達火山圖和差異基因Venn圖

通過對小地老虎轉錄組KEGG通路富集分析，結合Interperscan結構域預測和其他注釋信息得到小地老虎翅發育機制相關通路基因，并進行表達量分析與表達量差異比較分析，對差異表達基因繪制小地老虎Wnt信號通路基因表達熱圖（圖4）；并根據KEGG通路差異基因富集分析繪制小地老虎Wnt信號通路基因圖（圖5）。由小地老虎Wnt相關基因KEGG通路基因富集分析和表達量差異比較分析可知，Notum、Frizzled、CaN等53個元件均被注釋出，但是PRP、GBp、ICAT等17個元件沒有在小地老虎轉錄組中發現（圖4-5）。Wnt信號通路包含經典Wnt通路，非經典Wnt /鈣通路以及非經典planar cell polarity通路，其中非經典Wnt /鈣通路以及非經典planar cell polarity通路總基因數目較少，但是在小地老虎轉錄組中注釋出的基因較多，非經典Wnt/鈣通路只有一個元件未注釋出，非經典planar cell polarity通路只有兩個元件未注釋出；非經典Wnt /鈣通路，調控細胞內的鈣，非經典planar cell polarity通路調控細胞骨架，骨架決定細胞的形狀。

圖4 小地老虎Wnt信號通路基因表達熱圖

小地老虎胸部組織與頭部組織相比，Frizzled（Fz）、BMP and activin membrane-bound inhibitor（BAMBI）、（CKIα）、casein kinase 1（CKIε）、Axin、Myc proto-oncogene protein（c-myc）、E3 ubiquitinprotein ligase（Siah-1）、glypican 4（Knypek）、serine/threonine-protein phosphatase 2B catalytic subunit（CaN）、cullin 1（Cul1）等在小地老虎胸部組織中都上調表達，Notum、calcyclin binding protein（SIP）、classical protein kinase C alpha type（PKC）在胸部組織中下調表達。小地老虎胸部組織與腹部組織相比，Notum、Dally、Siah-1、CaN等在小地老虎胸部組織中都上調表達，wingless-type MMTV integration site family、Frizzled、Rac等在小地老虎胸部組織中下調表達。小地老虎胸部組織與小地老虎幼蟲組織相比，Frizzled、Dally、SIP、Cul1、Siah-1等上調表達，Rac、phosphatidylinositol phospholipase C beta（PLC）、Frizzled等下調表達（圖5）。

小地老虎胸部組織相對于頭部、腹部及幼蟲組織，只有Siah-1和CaN全部上調表達，Notum雖然在小地老虎胸部組織相對于頭部、幼蟲對比中上調表達，但在小地老虎胸部組織相對腹部組織中下調表達。SIP同時在小地老虎胸部相對于頭部和腹部比較中下調表達。

非經典planar cell polarity通路中，小地老虎胸部相對于頭部有表達差異的基因全部為表達上調，小地老虎胸部相對于腹部有表達差異的基因全部為表達下調。非經典planar cell polarity通路中，PLC在小地老虎胸部相對于頭部和幼蟲的基因表達中均為下調，下游元件CaN表現為表達上調。

圖5 小地老虎Wnt信號通路基因圖

3 討 論

Wnt信號通路相關基因在小地老虎幼蟲、小地老虎成蟲頭、胸、腹部均有表達。Wnt信號通路包含3種Wnt信號傳導通路：經典Wnt通路，非經典Wnt /鈣通路以及非經典planar cell polarity通路。非經典planar cell polarity通路是通過Wnt與Fz及其共受體結合而激活的。然后受體結合Dishevelled（Dsh），Dsh利用其PDZ和DIX結構域與morphogenesis 1（Daam1）的Dishevelled-associated激活因子1形成復合物。然后Daam1通過鳥嘌呤交換因子激活G-蛋白Rho［24］。Rho激活Rho相關激酶Rho-associated kinase（ROCK），這是細胞骨架的主要調節因子之一。Dsh還與Ras-related C3 botulinum toxin substrate 1（Rac1）形成復合物，并介導profilin與肌動蛋白的結合［25］。Rac1激活JNK，也可導致肌動蛋白聚合。丙氨酸與肌動蛋白結合可導致細胞骨架和原腸胚的重組。經典Wnt通路相關基因有一些在小地老虎轉錄組中未被注釋到，這可能與小地老虎的發育有關，有些基因沒有表達，也有可能與物種有關，小地老虎可能不含有某些基因。同時這些在小地老虎腹部表達的基因可能參與一些基因轉錄的調控，這也是此信號通路的主要功能。

Wnt信號通路相關基因在小地老虎幼蟲和成蟲不同組織表達差異揭示小地老虎可能的翅發育機制，非經典planar cell polarity通路調控細胞骨架，決定骨架細胞的形狀，Rac在小地老虎腹部相對小地老虎胸部上調表達，Rac也能間接影響JNK，從而影響一些基因的轉錄。非經典Wnt /鈣通路中，小地老虎胸部相對于頭部、腹部和小地老虎幼蟲，CaN、Siah-1都表現為上調表達，CaN在肌肉形成和功能的基本過程中起著重要的作用［26］，因此小地老虎成蟲胸部CaN的表達可能與其肌肉功能有關。本次試驗的小地老虎成蟲為處于遷飛狀態的小地老虎成蟲，胸部肌肉相關基因表達較高與其生物性狀相適應。

Wnt信號通路相關基因在小地老虎幼蟲、成蟲頭、胸、腹部的表達存在差異，小地老虎胸部組織與頭部組織、腹部組織和小地老虎幼蟲組織相比，CaN、Siah-1相關基因都顯著上調表達。這些基因有可能在其他昆蟲中也會有類似現象。對于昆蟲翅發育相關基因的研究主要集中于昆蟲羽化之前，本試驗表明，在昆蟲羽化后一些翅發育相關基因會繼續在成蟲胸部表達，也可能與翅部肌肉等組織的發育和代謝有關。這些對于昆蟲或者對于遷飛昆蟲的研究可能提供一些研究思路。近年來，隨著illumina、pacbio、nonapore等高通量測序技術的迅猛發展，以及CRISPR-Cas9技術的不斷發展，越來越多的科學問題通過基因組測序成功解決，也為研究昆蟲遷飛相關基因功能解析提供了新的技術方法［27］。

另外，影響昆蟲翅發育相關通路還有Decapentaplegic（Dpp）、Hedgehog、Hippo信號通路等。Dpp是參與果蠅發育的關鍵形態發生素，Dpp也是第一個經過驗證的分泌形態發生素［12-13］。Hedgehog（Hh）信號通路在控制細胞增殖、組織模式形成、干細胞維持和發育方面發揮著多種作用［15-18］。Hippo信號通路也被稱為Salvador或者Warts通路，通過調控細胞增殖和凋亡來控制動物的器官大小。該通路的名字來源于它的一個關鍵信號元件——蛋白激酶Hippo（Hpo）［19-20］。該基因的突變導致組織過度生長，或類似“河馬”表型。

昆蟲翅的進化對昆蟲的意義重大，探究影響昆蟲翅生長發育的相關通路，有助于了解昆蟲翅的進化及發育的過程。在眾多的昆蟲中，有相當一部分為農業害蟲，有些如棉鈴蟲、小地老虎等還是可進行遷飛的害蟲，探究昆蟲翅的發育機制，可以為揭開昆蟲遷飛面紗并研發新型防治農業害蟲方法提供幫助。本研究通過對小地老虎轉錄組的組裝及注釋，結合影響果蠅翅發育相關通路基因，找到許多潛在與小地老虎翅發育相關通路基因，為后續研究小地老虎翅發育相關基因提供了良好的基礎。

4 結論

通過對小地老虎轉錄組的組裝及注釋，最終組裝得到轉錄組243.4 Mb，組裝出轉錄本429 288個，獨立基因37 744個。小地老虎Wnt相關通路中，Notum、Frizzled、CaN等53個相關元件被注釋出。Wnt信號通路中的非經典planar cell polarity通路中，PlC在小地老虎胸部相對于頭部和幼蟲的基因表達中均為下調，下游元件CaN、Siah-1在小地老虎胸部相對于頭部、腹部及幼蟲組織的基因表達中均表現為表達上調，可能與其肌肉功能有關。