基于RNA-seq 技術的羅氏沼蝦(Macrobrachium rosenbergii)不同組織轉錄組比較分析*

李喜蓮 顧志敏 慎佩晶 徐 洋 張宇飛 高 強 程海華 陳雪峰

(農業農村部淡水漁業健康養殖重點實驗室 浙江省淡水水產遺傳育種重點實驗室 浙江省淡水水產研究所 湖州 313001)

羅氏沼蝦(Macrobrachium rosenbergii), 又稱淡水長臂大蝦、馬來西亞大蝦, 動物學分類上屬于節肢動物門(Arthropoda), 軟甲綱(Malacostraca), 十足目(Decapoda), 長臂蝦科(Palaemonidae), 沼蝦屬(Macrobrachium), 在各種類型的淡水和半咸水中都能生活, 易養殖且體型肥大, 肉質鮮美。因其適應能力強、生長速度快, 羅氏沼蝦已成為我國重要的淡水養殖蝦類之一。

近年來對羅氏沼蝦的研究主要包括遺傳多樣性研究, 性腺發育、基因克隆及生長等方面的研究。其中在多樣性研究中, 以微衛星為研究對象的如呂敏等(2019)研究異型雄性羅氏沼蝦遺傳多樣性, 戴習林等(2017)對種群 SSR 分析中樣本量及標記量對遺傳多樣性指標的影響進行分析; 或對各不同群體進行遺傳多樣性研究(董丁健等, 2020; 馮藝, 2018; 孫成飛等, 2015); 周曉敏等(2020)選取60 個SNP 位點對養殖群體和選育群體進行多樣性研究。羅氏沼蝦雌激素相關受體(Estrogen-related receptor, ERR)成為研究卵巢發育的一個重要基因(趙苗鑫等, 2017; 劉金磊,2018; 劉金磊等, 2018)。藥物對性腺分化及發育的影響研究主要包括壬基酚(薛海波, 2010; 李郁嬌, 2011;朱春華等, 2017; Guo et al, 2019)、十氯酮(Lafontaine et al, 2016, 2017)及三丁基錫(薛海波, 2010; 李郁嬌,2011; 吳維福等, 2013)這三種藥物, 以期研究對羅氏沼蝦卵黃蛋白原基因表達以及性腺發育的影響。Stalin 等(2019)研究了鈷-60 射線對淡水對蝦生殖障礙的影響; Tan 等(2019)研究了羅氏沼蝦性別逆轉與雄激素腺(AG)的關系。羅氏沼蝦基因研究集中在免疫相關(劉偉利等, 2017; 江婷佳, 2017; 邱慶慶等, 2019)、性別相關(俞炎琴, 2013; 姜建萍等, 2019; Abayed et al, 2019)、生長相關(葉成凱等, 2019; 邱慶慶, 2019;楊光等, 2020; Dong et al, 2020)和酶類(田榮等, 2016;張夏青等, 2016; 盧志杰等, 2019)。羅氏沼蝦生長研究集中在藥物對生長的影響(程安達等, 2019; 潘忠超等,2019; Tadese et al, 2020)、微生態制劑(朱光來等, 2019;趙臣澤, 2019)、生長環境(戴習林等, 2016; 陳建酬等,2017; 何竺柳, 2018; 朱其建等, 2019; 張俊功, 2019;Chen et al, 2019; Dong et al, 2020; Manickam et al,2020)、飼料中不同配比(楊樹浩等, 2018; 楊景豐等,2019; 張劍偉等, 2019; 單凡等, 2019; 黃黎明等,2019; Feng et al, 2019)和不同品系生長對比試驗(蔣飛等, 2013, 2014, 2016)。

隨著高通量測序技術的發展, 轉錄組技術也應用于羅氏沼蝦的研究中。郭梁等(2018)利用高通量測序技術和數字基因表譜等技術對感染螺原體的羅氏沼蝦進行免疫通路相關基因及其差異表達分析, 獲得轉錄本43405 個; 嚴賽峰等(2018)、李俊杰等(2018)和鄧澤森等(2018)對感染螺原體的羅氏沼蝦高通量測序結果開展了 SSR 位點、SNP 位點和重要免疫通路相關 microRNA 研究和分析; 李威霖(2018)對羅氏沼蝦肢體再生組織進行了轉錄組及其差異表達基因研究, 得到 Unigenes 總數目為 87783; 王傳聰等(2018)對羅氏沼蝦肝胰腺組織轉錄組測序結果進行 SSR 檢測與分析, 獲得15356 個SSR 位點; 陳雪峰等(2019)采用Illumina HiSeqTM4000 高通量測序研究羅氏沼蝦卵巢發育四個時期卵巢組織的差異, 卵巢發育四個時期共獲得 95379 個 Unigenes。Jiang 等(2019a)通過高通量測序技術研究了羅氏沼蝦雌性(ZW)、雄性(ZZ)和超雌個體(WW)的性腺差異。Pasookhush 等(2019)采用轉錄組技術分析了羅氏沼蝦幼蝦感染新型冠狀病毒的反應。Jiang 等(2019b)對未成年雌蝦和成年雌蝦的眼柄組織進行轉錄組分析, 共獲得 53878 個Unigenes; Cao 等(2017)研究羅氏沼蝦正常樣本和WSSV 感染樣本的淋巴組織轉錄組, 分別獲得73658和72374 個Unigenes; Rao 等(2016)分別研究了羅氏沼蝦正常樣本和 WSSV 感染樣本的肝胰腺轉錄組, 共獲得 63584 個 Unigenes。

本研究以羅氏沼蝦7 個組織(眼柄、肝臟、卵巢、鰓、心臟、肌肉、精巢)為研究對象, 采用新一代高通量測序技術進行轉錄組測序分析, 通過 Trinity 軟件組裝、數據庫功能注釋、基因表達差異分析和SSR、SNP 位點篩選, 獲得羅氏沼蝦分子遺傳信息, 以期為進一步研究羅氏沼蝦遺傳多樣性、功能基因及基因表達差異提供理論數據, 同時為深入研究羅氏沼蝦生長阻滯現象提供一定的數據基礎。

1 材料與方法

1.1 實驗材料

實驗用羅氏沼蝦(Macrobrachium rosenbergii)樣品取自浙江省淡水水產研究所八里店綜合試驗基地,體重為(5.34±1.22) g, 體長為(6.13±0.42) cm。樣品經解剖后取各組織(眼柄、肝臟、卵巢、肝臟、心臟、肌肉、精巢), 共7 個組織。每個組織取自3 個個體(見表 1)。將分裝好的樣品迅速投入液氮中速凍, 置于-80°C 冰箱中保存, 干冰運輸。

表1 組織及樣品名稱列表Tab.1 List of tissue and sample names

1.2 方法

1.2.1 RNA 提取 各組織樣品均選用3 個個體的混合樣組織, 使用 TRIzol?試劑按照制造商的說明從各組織中提取總 RNA (Invitrogen), 使用 DNase I(TaKaRa)去除基因組DNA。使用Nanodrop ND-2000分光光度計(美國賽默飛)、Aglient 2100 分析儀器對總RNA 的純度、濃度和完整性進行檢測。RIN 值＞7的RNA 用于下游實驗。

1.2.2 mRNA-seq 文庫構建和Illumina 測序 使用mRNA-seq 樣品制備試劑盒(Illumina, San Diego, CA)按照試劑說明書步驟構建了 mRNA-seq 文庫和Illumina 測序文庫。

1.2.3 質量控制和從頭轉錄組組裝 使用Fast QC程 序 (http://www.bioinformatics. babraham.ac.uk/projects/fastqc/)檢查Illumina Hiseq 2500 平臺在校準前產生的原始序列讀取的質量。低質量讀數低于閾值質量20; 將長度小于50 bp 的reads 以及包含適配序列、ploy-N 和來自原始數據的測序引物的reads 去除,得到干凈的reads同時, 對質控數據的誤差率%、Q30、GC-含量%和sequence 重復水平進行了評價。所有的后續分析和注釋都依賴于高質量的clean reads。

使用Trinity (v 2.8.5)軟件(https://sourceforge.net/projects/trinityrnaseq/)對質控后的數據進行從頭組裝,獲得轉錄本序列和Unigenes 序列。最后, 只有長度≥300 bp 的unique contigs 才能用于組裝后的下游研究。

1.2.4 生物信息學分析 功能注釋。拼接得到的Unigenes 序列, 使用 BLASTX 比對(BLAST+2.7.1, 比對標準: E 值不大于 1e-5)與非冗余蛋白數據庫(Non-Redundant Protein Sequence Database, NR)、SWISS-PROT 蛋白質序列數據庫(Manually Annotated and Reviewed Protein Sequence Database, SWISS-PROT)、基因本體數據庫(Gene Ontology, GO)、直系同源蛋白數據庫(Cluster of Orthologous Groups of Proteins,COG), 京都基因和基因組百科全書 KEGG 數據庫比對。使用Trinity 軟件自帶的ORF 預測模塊進行開發閱讀框(ORF)預測。

1.2.5 差異表達分析 通過 Bowtie2 軟件將 reads映射到組裝的 Unigenes 序列, 使用 RSEM 軟件根據比對結果來計算特定轉錄本的表達量水平。衡量基因表達水平的標準為 RPKM 值(Reads per kilobase of exon model per million mapped reads), 即每一百萬條序列中, 每個基因以一千個堿基為單位, 比對上的reads 個數。

使用DE-Seq 軟件分析各個組織之間的差異表達,從而找到差異基因組。以差異倍數(Foldchange)＞ 2 和假發現率(FDR)調整顯著性值≤0.05 為判斷Unigenes表達顯著性的閾值。利用GO、egg NOG、KEGG 或thology (KO)和KEGG 通路富集分析對DEGs 進行分類, 并對潛在的生物學途徑中的DEGs 進行評價。GO和KEGG 途徑富集分析中P 值小于0.05 的過程、功能或成分在DEGs 中被認為存在顯著差異。

2 結果

2.1 原始測序數據質控及統計

羅氏沼蝦7 個不同組織測序共得到344151324 條原始序列, 質控后得到有效RNA-seq 311475706 條。每個個體的測序量為6.65—9.04 Gb, 平均7.38 Gb。去除接頭序列, 截去連續4 個堿基平均質量值低于20的部分, 舍去長度少于 50 bp 的 reads。共得到311475706 條reads。過濾后每個個體的測序量為5.3—8.45 Gb, 平均6.60 Gb(見表2)。每個個體的測序量為每4 個堿基的平均質量值均大于20。

表2 質控后數據統計表Tab.2 Statistics of data after quality control

2.2 轉錄本拼接

轉錄本拼接后共獲得 95220 個 Unigenes, 總Unigenes 長度為 101401098 bp。總的轉錄本數量為145717 個, 總的轉錄本長度為 207379988 bp。每個Unigenes 的平均長度為 1064.9 bp, 最長的 Unigenes長為36137 bp, N50 值為1553。大部分功能基因長度分布在401—600 bp, 占比達到34.27%(見圖1)。

2.3 功能注釋

在NR、GO、COG、KEGG、SWSS-PROT 五個公共數據庫中對獲得的 95220 個 Unigenes 進行功能注釋, 其余未在NR 數據庫中找到的Unigenes 可能為新的蛋白。

根據GO 數據庫, 總共有18485 個基因被歸類到三個主要的GO 類別中: 生物過程、分子功能和細胞成分。其中, 以“細胞”(14938)、“細胞部分”(14925)、“細胞器”(10827)、“生物調節”(9861)為主(見圖 2)。

根據COG 功能分類分成26 類, 其中數量前三位的為功能未知 2172 個(13.56%), 信號轉導機制 2070個(12.92%), 翻譯后修飾, 蛋白質轉換, 伴侶1660 個(10.36%)(見圖 3)。

圖1 組裝序列長度分布圖Fig.1 Assembly sequence length distribution diagram

圖2 GO 統計二級統計圖Fig.2 Secondary chart of GO level 2

五個數據庫 NR、GO、COG、KEGG、SWSS-PROT分別注釋到 19881、18485、15798、9147、13684 個Unigenes, 在 NR 數據庫中比對注釋的信息最多, 注釋19881 個Unigenes, 占比達20.88%。在NR、GO、COG、KEGG、SWSS-PROT 五個數據庫中都注釋到的 Unigenes 有 7848 個(見圖 4)。

將基因根據參與的 KEGG 代謝通路分為 5 個分支(見圖5): 細胞過程(A, Cellular Processes), 環境信息處理(B, Environmental Information Processing), 遺傳信息處理(C, Genetic Information Processing), 代謝(D, Metabolism), 有機系統(E, Organismal Systems)。以上5 個分支中含量的最多的類型分別為: 全局和概率地圖(Global and overview maps), 轉化(Translation),單組織過程(Signal transduction), 運輸和分解代謝(Transport and catabolism)、內分泌系統(Endocrine system)。

2.4 差異基因表達分析

圖3 COG 分類統計圖Fig.3 COG classification statistical chart

圖4 注釋信息統計韋恩圖Fig.4 Annotate the information with statistical Venn diagrams

7 個不同組織鑒定到的mRNA 數量由多到少排序為: G(鰓) ＞ T(精巢) ＞ H(心臟) ＞ L(肝臟) ＞ E(眼柄)＞ O(卵巢) ＞ M(肌肉)(見圖 6)。在 7 個組織中共同表達的基因數為15260 個。

2.5 KEGG 通路分析

對羅氏沼蝦轉錄組的 KEGG 分析顯示, 9148 個Unigenes 被注釋到330 條KEGG 通路中, 其中代謝途徑(1578 個)、核糖體(540 個)、次生代謝物的生物合成(528 個)通路數量居前三(見圖7)。

330 條通路中, 其中信號通路包括: PI3K-Akt 信號通路、Rap1 信號通路、催產素信號通路、Hippo信號通路、MAPK 信號通路、胰島素信號通路、雌激素信號通路、Ras 信號通路、甲狀腺激素信號通路、鈣信號通路、HIF-1 信號通路、FoxO 信號通路、生成信號通路、Wnt 信號通路、PPAR 信號通路、促性腺激素信號通路、趨化因子信號通路、Adipocytokine信號通路、點樣受體信號通路、ErbB 信號通路、mTOR信號通路、Jak-STAT 信號通路、T 細胞受體信號通路、p53 信號通路、腫瘤壞死因子信號通路、VEGF 信號通路、B 細胞受體信號通路、Notch 信號通路、催乳激素信號通路、NF-kappa B 信號通路、Fc epsilon RI信號通路等。

2.6 SSR 分析

在得到的95220 條序列中篩選SSR 位點共找到SSR位點37751 個, 這些位點存在于25520 條序列中, SSR 發生頻率為 26.80%。其中單條序列中含多個SSR 位點的序列有7962 條, 復合型的SSR 位點3384 個。SSR 位點中單堿基重復、二堿基重復、三堿基重復、四堿基重復、五堿基重復、六堿基重復個數分別為 14919(39.52%)、14715(38.98%)、7577、488、32 和 20, 其中單堿基和二堿基重復含量居第一和第二位(見圖8)。

2.7 SNP 分析

對獲得的序列進行 SNP 位點分析, 共發現3228575 個 SNP 標記(見表 3), 其中包括 C:G-＞A:T、C:G-＞G:C、C:G-＞T:A、T:A-＞A:T、T:A-＞C:G、T:A-＞G:C6中類型, 其中T:A-＞C:G 和C:G-＞T:A 這兩種堿基替換占比較高, 分別占總數的 33.45%和 33.36%。7 個不同組織中發現的SNP位點數量在鰓(23.12%)上最高, 其次是心臟(17.89%)、精巢(17.51%)和肝臟(16.28%), 在肌肉組織中的SNP 位點數量最少(5.43%)。

圖 5 KEGG 注釋Fig.5 The KEGG annotation

圖6 每個樣本鑒定得到的mRNA Upset 圖Fig.6 The upstate map of mRNA from each sample

圖 7 KEGG 通路列表(前20)Fig.7 The KEGG pathway list (Top 20)

圖8 SSR 位點的分布情況Fig.8 The distribution of SSR sites

3 討論

3.1 轉錄組技術在其他蝦上的應用

隨著轉錄組技術的迅速發展, 轉錄組技術被應用在水產動物的研究中。其中蝦類轉錄組研究主要集中在差異基因的篩選和候選基因的發掘上。

紅、黃和透明3 種純色米蝦皮膚組織的轉錄組獲得 45434 條 Unigenes(林師, 2017)。波紋龍蝦肝胰腺和卵巢組織進行了轉錄組測序獲得Unigenes 共74124個(李斌等, 2016)。脊尾白蝦成蝦樣品及胚胎樣品共獲得 47574 條 Unigenes(孫政, 2012)。

表3 不同SNP 類型在各組織之間的分布情況Tab.3 Distribution of different SNP types among different tissues

紅螯螯蝦肝臟、精巢和卵巢組織共獲得了67369個 Unigenes(李喜蓮等, 2019); 次級卵黃發生時期卵巢組織共得到69261 條Unigenes(康鵬飛, 2017)。

日本沼蝦正常性成熟的和性早熟的卵巢組織中共獲得63336 個Unigenes(江紅霞, 2017); 感染白斑綜合征病毒(WSSV)個體的肝胰腺轉錄組共獲得 64049個 Unigenes(趙才源等, 2018); 亞硝酸鹽脅迫下肝臟共獲得19022 個Unigenes, 氨氮脅迫下肝臟組織共獲得 63453 個 Unigenes (于杰倫, 2019)。

凡納濱對蝦(Litopenaeus vannamei)五個早期發育時期(卵裂期、囊胚期、原腸期、肢芽幼體期發育至膜內無節幼體期)轉錄組測序, 共得到 66815 條Unigenes(隗健凱, 2015); 性腺組織獲得 Unigenes 65218 個(韋嬪媛, 2017); 低溫脅迫下肝胰腺組織獲得50921 條 Unigenes (董麗君等, 2019)。

本研究從正常個體 7 個組織轉錄組文庫中共獲得99520 個Unigenes, 這個數量較前人對于蝦類轉錄組研究得到的 Unigenes 數量都高, 這可能與分析的組織數量大、覆蓋面廣有一定的關系。

3.2 功能注釋

本研究從七個組織共獲得 Unigenes 99520 個,N501553, 平均長度為1064.9 bp。在五個數據庫中注釋到的Unigenes 共有20368 個, 占到總數的21.39%,這比前人研究的轉錄組結果Unigenes 的注釋率都低。肝胰腺注釋率31%和35.31%, 卵巢31.46%和54.44%,眼柄29.3%, 淋巴器官29.46%, 再生肢體基部注釋率37.23%。估計與本研究中首次對羅氏沼蝦的心臟、肌肉、精巢和鰓等組織進行轉錄組高通量測序, 引入了較多在這幾個組織中特異表達的基因, 這些新獲得的序列在以上五個數據庫中得不到注釋, 從而降低了注釋率。還有部分序列未能被注釋, 可能與序列的長度有關, 過短的序列也會造成無法注釋和分類; 無法注釋的另一種可能是近緣物種序列信息的缺乏,導致無法通過同源序列比對得到注釋。

3.3 KEGG 通路

對羅氏沼蝦轉錄組的 KEGG 分析顯示, 9148 個Unigeness 被注釋到 KEGG 數據庫中并分布在330 個已知途徑中, 與免疫相關的通路如FoxO 信號通路和Rap1 信號通路。FoxO 轉錄因子家族在細胞生理事件中調控基因的表達, 包括凋亡、細胞周期控制、葡萄糖代謝、氧化應激抵抗和壽命。FoxO 蛋白的一個主要調控機制是對磷脂酰肌醇3 激酶(PI3K)下游的絲氨酸蘇氨酸激酶Akt/蛋白激酶B(Akt/PKB)進行磷酸化,這是對胰島素或幾種生長因子的反應。FoxO 信號通路在本研究結果中涉及 93 個 Unigeness。Rap1 是一種小型GTPase, 它控制多種過程, 如細胞黏附, 細胞-細胞連接的形成和細胞極性。Rap1 通過調節多種細胞類型中整合素等黏附分子的功能, 在細胞-細胞和細胞-基質相互作用的控制中發揮主導作用。Rap1 還以高度依賴于細胞類型的方式調控 MAPK 活性。Rap1 信號通路在本研究中共涉及201 個Unigeness。這些結果的獲得都將為進一步研究羅氏沼蝦抗性相關基因提供科學依據。

4 結論

本研究對羅氏沼蝦 7 個不同組織的轉錄組數據進行分析, 獲得99520 個Unigenes, 這比單一組織或較少組織受到環境脅迫或細菌、病毒感染獲得數據具有較高的可信度, 這一結果將大大豐富羅氏沼蝦的基因數據庫資源。與此同時, 20368 個Unigenes 在五大數據庫中得到注釋。各個組織間差異基因也得到進一步的分析, 本研究還篩選得到大量的 SSR 位點和SNP 位點, 這些分子標記也將在后續羅氏沼蝦分析標記的開發中起到重要的作用, 為進一步深入挖掘和開發利用羅氏沼蝦功能基因提供參考。