利用OxfordNanopore長讀長RNA-Seq對發(fā)育中橄欖果蠅胚胎的轉錄組學研究

英文題目：Transcriptome landscape of the developing olive fruit fly embryo delineated by Oxford Nanopore long-read RNA-Seq

中文題目：利用OxfordNanopore長讀長RNA-Seq對發(fā)育中橄欖果蠅胚胎的轉錄組學研究

研究背景

橄欖果蠅或橄欖蠅（Bactroceraoleae）是栽培橄欖樹主要的害蟲，像所有的昆蟲一樣，橄欖蠅會完全變態(tài)。然而，尚未探索早期胚胎發(fā)育過程中發(fā)生的轉錄動力學，而在沒有完全注釋的基因組的情況下進行詳細的轉錄組學分析具有挑戰(zhàn)性。收集發(fā)育前6個小時每小時的橄欖蠅胚胎進行ONT測序，每個時間點獲得3100萬reads，與橄欖蠅基因組比對效率在98%，全長覆蓋率大于50%，在發(fā)育的前六個小時中檢測到68％的預測的基因的表達。鑒定了3553個新基因和共79,810個轉錄本，與NCBI預測的轉錄組相比，轉錄組多樣性增加了四倍。胚胎發(fā)育的前六個小時的特征在于顯著的轉錄組變化，每個胚胎的轉錄物總數從胚胎發(fā)育的第一小時到第二小時降至一半。基于時間共表達的基因聚類，在胚胎發(fā)育的前六個小時表達的基因的基因集富集分析顯示參與轉錄和翻譯，大分子生物合成和神經發(fā)育的基因高度富集?？傊?，cDNA分子的全長測序詳細表征同種型復雜性和B.oleae的第一胚胎階段的轉錄動力學。

結果分析

1、目前B.oleae基因組組裝和基因組注釋

橄欖蠅有六對染色體，其中包括一對異色性染色體，雄性為異性性染色體，最初通過qPCR估計B.oleae基因組大小在322Mb左右。作者之前提交過一版本基因組，該注釋包含總共13936個基因和假基因，其中,分別被預測為蛋白質編碼13198，非編碼392和假基因346。此外，預測2,759個基因具有可變剪接。而基因和轉錄本的平均長度分別為9,597bp和2,259bp，長的基因為497,921bp，而長的轉錄本為59,475bp。

2、橄欖蠅的轉錄組測序鑒定出新的基因和亞型

作者獲得了3100萬次reads，其中使用Canu對2200萬reads（71％）進行了錯誤糾正，只關注全長reads，通過確定為具有5’引物和poly（A）和3’引物的序列全長序列。ToFU轉錄組結果包含總共11,883個基因和79,810個亞型（表北京百邁客生物科技有限公司1），其中8330個基因與NCBI注釋基因匹配，而3553個基因是新的。所有這些都對應于橄欖蠅轉錄組在同種型水平上比當前NCBI注釋的四倍擴增。針對UniprotSwiss-prot數據庫的預測蛋白質序列的搜索顯示高比對的是雙翅目，其次是哺乳類動物（圖1B）。盡管鑒定了3553個新基因，但預計只有269個含有開放閱讀框。與注釋基因相比，超過50％的新基因是單外顯子，其中超過80％的基因是多外顯子并且還包含更高百分比的單同種型基因。且新基因的表達量要低于已注釋基因。在結構上，SQANTI根據其剪接點和基因組坐標將轉錄本分為9類，分別為FSM，ISM，NIC，NNC，基因組基因，反義，融合，基因間區(qū)和基因內含子（圖1E）。

3、RNA-Seq數據的直接標準化優(yōu)于相對標準化

為了獲得轉錄數，在cDNA合成步驟中以每個時間點使用的每個胚胎數的恒定比率添加ERCC內部加標RNA標準。在測序和比對后，通過2個步驟實現標準化，使用Mandalorion對測序深度進行相對標準化，得到每10000個映射讀數的每個基因讀數的轉錄本豐度（RP10K），我們注意到我們的RNA標準品的相對標準化豐度隨時間變化，這很可能是由于胚胎中poly（A）RNA含量的變化（圖2A）。在第二步中，我們使用ERCC標準生成的標準曲線將我們的相對計數轉換為每個胚胎的計數。在這里，我們注意到ERCC標準并沒有隨著時間點的不同而發(fā)生顯著的變化（圖2B）。有趣的是，與相對表達相反，當所有基因的每個胚胎的轉錄物的數量相加并在時間點上繪圖時（圖2C和D），該譜與每個胚胎產生的cDNA相似（補充圖4A），因此驗證標準化方法。相鄰時間點北京百邁客生物科技有限公司顯示出比遠距離時間點更高的基因表達相關性，其中連續(xù)樣本的Spearman相關性始終等于或高于0.96。

4、胚胎的總mRNA含量和生物學重復驗證絕對定量

通過將每個胚胎的每個基因的所有轉錄物相加并計算以毫微克計的當量來計算每個胚胎的總mRNA（圖3A），總mRNA在1hpo時從1.2ng/胚胎下降至2hpo時為0.61ng/胚胎，然后在3hpo時增加至1.49ng/胚胎，然后在6hpo采樣結束時降至0.93ng/胚胎，反映了每個胚胎的總轉錄本（圖2D）。假設2-5％的總RNA是多腺苷酸化的，mRNA水平與我們在每個胚胎中獲得的總RNA產量一致。

進一步尋求確定觀察到的表達模式是否可以在不同的生物樣品組中重復并使用實時定量PCR（qPCR），這是當前用于量化基因表達的標準方法。然而，qPCR基因表達的相對方法需要鑒定內參基因其表達在樣品間保持穩(wěn)定。使用作者的數據來確定這些常見內參基因以及GAPDH的表達水平的變化，評估了3個基因的qPCR表達;SRY，HID和LINGERER，用RPL19和14-3-3zeta作為內參基因，在不同的生物重復樣品組中（跳過6小時時間點）。我們觀察到類似的基因表達趨勢（特別是14-3-3zeta與RPL19相比）。為了進一步探索表達譜，我使用差異表達的明顯的基因進行主成分分析（PCA）和層次聚類分析，第一個主成分將前3個時間點與最后3個時間點分開，分層聚類進一步表明前3個和后3個時間點是分別共同聚類的。

5、母體到受精卵轉變表明胚胎mRNA含量發(fā)生了顯著變化

線蟲，昆蟲，魚類，兩棲動物和哺乳動物在內的許多后生動物的發(fā)育胚胎的特征在于顯著的轉錄變化，其中之一是胚胎依賴從母體到受精卵轉錄物的變化。MZT過程包括兩個階段，首先是在卵子發(fā)生過程中清除大部分母本轉錄本和最初加載到卵母細胞中的蛋白質，然后開始受精卵轉錄。在黑腹果蠅中，MZT已被廣泛研究，胚胎依賴于母體轉錄本和蛋白質直至受精后2-3小時。然而，在MZT期間，在2hpf結束時，母系編碼的蛋白質中有多達20％的母系供應的轉錄物不穩(wěn)定，而另外15％的母本轉錄本通過化學編碼的蛋白質3hpf不穩(wěn)定。黑腹果蠅母系不穩(wěn)定基因富含細胞周期功能，而母系穩(wěn)定基因則富含家族保持功能，如代謝，翻譯。

使用我們的時間進程數據來闡明MZT在B.oleae中的作用機制，這個過程在我們所知的范圍之前還沒有被研究過。我們在檢測跨時間點的發(fā)育過程中每個胚胎的總mRNA含量時發(fā)現了一個有趣的現象（圖3A），與1hpo時相比，2hpo時每胚胎的總mRNA下降51％，3hpo時相比2hpo時增加143％。實際上，表達基因的數量在1,2和3hpo之間是相似的，使用時間進程的數據，利用GFOLD在連續(xù)時間點之間進行差異表達，GFOLD是為沒有生物重復的樣品設計的，已顯示GFOLDlog2倍數變化與qPCR確定的倍數變化相關性較好。使用±0.5的Gfold截止值將基因編碼為上調或下調。確定了1496個基因，這些基因在2hpo時比1hpo下調。這些基因富含母體降解的轉錄物，在此稱為母體降解的基因。實際上，在1hpo時，母體降解基因的表達水平高于其他基因的表達水平，與其他基因相比，相同基因在2hpo時表現出相似的表達水平，表明這些基因不穩(wěn)定至其他基因的基礎水平。

我們對3個類別基因集進行富集分析;母體降解基因，受精卵基因和母體穩(wěn)定/上調基因。母體降解基因也是1hpo中表達高的基因，富含細胞過程，發(fā)育和新陳代謝（圖6）。母體降解的基因也富含轉錄因子，例如DREF，BEAF-32A，PNR，它們是相應的果蠅同源基因。類似母體降解基因，母體穩(wěn)定/上調基因在翻譯，生物合成過程，基因表達，代謝過程等方面得到豐富，反映了快速生長胚胎的高代謝活性，DREF轉錄因子也在這些基因中富集。受精卵基因富含特定的組織形成和發(fā)展過程，包括：后腸發(fā)育，模式規(guī)范，消化道形態(tài)發(fā)生等功能。

6、基于時間表達動態(tài)的基因聚類

基因表達是一個嚴格調控的過程，在胚胎發(fā)育過程中，基因表達的時空動態(tài)調控對器官的正常發(fā)育至關重要?；谄鋾r間表達動力學的基因聚類不僅將表達基質的復雜性降低為簡單的基因集，而且還可以鑒定具有與先前所示相似的生物學功能的基因。實際上，我們鑒定了在不同時間點表達達到峰值的基因，證明了高度動態(tài)的轉錄本動力學，并表明這些基因在確定的發(fā)育期間具有特定的作用（圖7）。我們進一步將這些集群分為3組;1）基因在3hpo達到峰值并且通常隨時間降低，稱為早期基因（圖7A），2）基因，其表達維持3-5hpo，稱為中間基因（圖7B和C），和3）基因其表達僅在5和/或6hpo時增加，稱為晚期基因（圖7D）。富集顯示，如先前在母體類別的基因中所觀察到的，早期基因和中間基因在細胞過程和代謝過程中富集，而晚期基因在專門的發(fā)育過程中富集。

7、長讀長RNA-Seq完善了基因在性別決定途徑中的注釋

雙翅目昆蟲的性別決定機制在很大程度上是相同的，在果蠅（Drosophilamelanogaster）中，性別決定機制已被廣泛研究，性致死基因（sxl）作為主要調節(jié)因子，根據性染色體與常染色體的比例，調節(jié)本身和變異基因（tra）的性別特異性選擇性剪接。tra反過來調節(jié)性別特異性雙重性別的選擇性剪接（dsx），這是級聯(lián)的最后一個成員和性別差異發(fā)育的中介。B.oleae同源物已經鑒定了sx1，tra和dsx。然而，主要的調控方式仍然難以捉摸。已經提出橄欖蠅的性別決定在胚胎發(fā)育的前6小時內發(fā)生，并且通過類似于黑腹果蠅中的性別決定機制的tra和dsx的可變剪接來調控。在轉錄組數據中，我們能夠觀察到tra和dsx的各種可變剪接。使用來自成年雄性和雌性頭部的數據來識別性別特異性亞型。在dsx的情況下，與成齡相比，早期發(fā)育階段的亞型復雜性顯著不同。在數據中，看到具有不同轉錄起始位點和長度較長的亞型，因為在發(fā)育的早期胚胎階段存在突出的亞型（圖8），成齡頭部組織中的這些亞型轉移到較短的亞型中，外顯子4存在于雌性中但不存在于雄性中。由于我們無法在早期胚胎階段檢測到這些性別特異性亞型，我們認為它們的表達在發(fā)育過程后期開始。然而，它們在后期階段的積累代表了在發(fā)育早期階段激活的性別調控系統(tǒng)。北京百邁客生物科技有限公司。

討論

由于受精卵的興趣和令人費解的性質，胚胎發(fā)育已在模型生物中廣泛研究。然而，在生物中沒有很好表征，長讀長RNA-seq有可能揭示迄今未知的基因，完善注釋信息，并擴展亞型多樣性。作者匯集了混合性別橄欖果蠅（Bactroceraoleae）胚胎，這些胚胎在產卵的前6個小時后每小時收集一次。這些時期的研究很有意義，因為之前對B.oleae的研究表明，在此期間開始通過基因的可變剪接介導的性別決定機制。此外，來自地中海實蠅（Ceratitiscapitata）的證據表明，在此期間發(fā)生極細胞（原始生殖細胞）的建立，因此，在此期間闡明轉錄狀態(tài)非常重要。

作者在轉錄組中包括雄性和雌性頭部以擴展轉錄組數據?？偟膩碚f，我們的轉錄分析工作提供了豐富的資源來識別早期發(fā)育基因和轉錄異構體以及一系列廣泛的可變剪接變體。