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

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

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

近年來對羅氏沼蝦的研究主要包括遺傳多樣性研究, 性腺發(fā)育、基因克隆及生長等方面的研究。其中在多樣性研究中, 以微衛(wèi)星為研究對象的如呂敏等(2019)研究異型雄性羅氏沼蝦遺傳多樣性, 戴習(xí)林等(2017)對種群 SSR 分析中樣本量及標(biāo)記量對遺傳多樣性指標(biāo)的影響進(jìn)行分析; 或?qū)Ω鞑煌后w進(jìn)行遺傳多樣性研究(董丁健等, 2020; 馮藝, 2018; 孫成飛等, 2015); 周曉敏等(2020)選取60 個SNP 位點對養(yǎng)殖群體和選育群體進(jìn)行多樣性研究。羅氏沼蝦雌激素相關(guān)受體(Estrogen-related receptor, ERR)成為研究卵巢發(fā)育的一個重要基因(趙苗鑫等, 2017; 劉金磊,2018; 劉金磊等, 2018)。藥物對性腺分化及發(fā)育的影響研究主要包括壬基酚(薛海波, 2010; 李郁嬌, 2011;朱春華等, 2017; Guo et al, 2019)、十氯酮(Lafontaine et al, 2016, 2017)及三丁基錫(薛海波, 2010; 李郁嬌,2011; 吳維福等, 2013)這三種藥物, 以期研究對羅氏沼蝦卵黃蛋白原基因表達(dá)以及性腺發(fā)育的影響。Stalin 等(2019)研究了鈷-60 射線對淡水對蝦生殖障礙的影響; Tan 等(2019)研究了羅氏沼蝦性別逆轉(zhuǎn)與雄激素腺(AG)的關(guān)系。羅氏沼蝦基因研究集中在免疫相關(guān)(劉偉利等, 2017; 江婷佳, 2017; 邱慶慶等, 2019)、性別相關(guān)(俞炎琴, 2013; 姜建萍等, 2019; Abayed et al, 2019)、生長相關(guān)(葉成凱等, 2019; 邱慶慶, 2019;楊光等, 2020; Dong et al, 2020)和酶類(田榮等, 2016;張夏青等, 2016; 盧志杰等, 2019)。羅氏沼蝦生長研究集中在藥物對生長的影響(程安達(dá)等, 2019; 潘忠超等,2019; Tadese et al, 2020)、微生態(tài)制劑(朱光來等, 2019;趙臣澤, 2019)、生長環(huán)境(戴習(xí)林等, 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)。

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

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

1 材料與方法

1.1 實驗材料

實驗用羅氏沼蝦(Macrobrachium rosenbergii)樣品取自浙江省淡水水產(chǎn)研究所八里店綜合試驗基地,體重為(5.34±1.22) g, 體長為(6.13±0.42) cm。樣品經(jīng)解剖后取各組織(眼柄、肝臟、卵巢、肝臟、心臟、肌肉、精巢), 共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 的純度、濃度和完整性進(jìn)行檢測。RIN 值＞7的RNA 用于下游實驗。

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

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

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

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

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

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

2 結(jié)果

2.1 原始測序數(shù)據(jù)質(zhì)控及統(tǒng)計

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

表2 質(zhì)控后數(shù)據(jù)統(tǒng)計表Tab.2 Statistics of data after quality control

2.2 轉(zhuǎn)錄本拼接

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

2.3 功能注釋

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

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

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

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

圖2 GO 統(tǒng)計二級統(tǒng)計圖Fig.2 Secondary chart of GO level 2

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

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

2.4 差異基因表達(dá)分析

圖3 COG 分類統(tǒng)計圖Fig.3 COG classification statistical chart

圖4 注釋信息統(tǒng)計韋恩圖Fig.4 Annotate the information with statistical Venn diagrams

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

2.5 KEGG 通路分析

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

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

2.6 SSR 分析

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

2.7 SNP 分析

對獲得的序列進(jìn)行 SNP 位點分析, 共發(fā)現(xiàn)3228575 個 SNP 標(biāo)記(見表 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 這兩種堿基替換占比較高, 分別占總數(shù)的 33.45%和 33.36%。7 個不同組織中發(fā)現(xiàn)的SNP位點數(shù)量在鰓(23.12%)上最高, 其次是心臟(17.89%)、精巢(17.51%)和肝臟(16.28%), 在肌肉組織中的SNP 位點數(shù)量最少(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 轉(zhuǎn)錄組技術(shù)在其他蝦上的應(yīng)用

隨著轉(zhuǎn)錄組技術(shù)的迅速發(fā)展, 轉(zhuǎn)錄組技術(shù)被應(yīng)用在水產(chǎn)動物的研究中。其中蝦類轉(zhuǎn)錄組研究主要集中在差異基因的篩選和候選基因的發(fā)掘上。

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

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

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

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

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

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

3.2 功能注釋

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

3.3 KEGG 通路

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

4 結(jié)論

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