基于單細胞轉錄組測序技術篩選蟾蜍耳后腺免疫修復基因

關鍵詞：單細胞轉錄組測序；中華大蟾蜍；耳后腺；免疫修復基因；實時熒光定量PCR

中圖分類號：S865.9：Q781 文獻標識號：A 文章編號：1001-4942（2024）11-0133-08

蟾酥是蟾蜍科動物中華大蟾蜍（Bufo bufo gargarizans Cantor）或黑眶蟾蜍（Bufo melanostictusSchneider）耳后腺和皮膚腺體的干燥分泌物，具有抗癌、強心、鎮痛、抗炎、抑菌等多種藥理作用。蟾酥是多種臨床藥物，如牙痛一粒丸、六神丸、麝香保心丸等的原材料，每年市場需求量極大，且價格逐年增長。隨著野生中華大蟾蜍種群資源的減少，蟾酥藥材需求量增大，人工馴養繁殖中華大蟾蜍獲取蟾酥藥材，成為行業熱點。然而，人工飼養蟾蜍取酥后，因耳后腺損傷感染導致蟾蜍死亡率極高，大多蟾蜍在進行一次采酥后就會死亡，這對于蟾蜍資源來說是一種極大的浪費，且蟾酥產量也達不到市場需求。

高通量測序技術已被廣泛用于農業研究，包括動植物全基因組測序、基因組重測序、轉錄組測序等。其中轉錄組測序技術通過與代謝組學等技術結合已被用于多種植物的機制研究。在轉錄組測序技術的基礎上，單細胞轉錄組測序（single-cell RNA-sequencing，scRNA-seq）技術是可以在單個細胞水平上對RNA進行高通量測序和分析的一項新技術。該技術可以同時分析成千上萬個細胞，能很好地將群體內單個細胞異質性區分開來，顯著提高了我們對動物組織、器官和有機體復雜性的理解。單細胞測序技術已應用于多種動物研究，張衛東等在絨山羊皮膚組織中鑒定出17個主要細胞類群，篩選出3個在絨山羊毛乳頭細胞中特異性表達的基因SOX2、FGF7與APOD;崔寶山等對牛體內雌性和雄性囊胚進行單細胞測序研究，發現二者的轉錄組信息存在巨大差異，并指出影響胚胎發育性別的特異性基因為FBP1、GADD45G、FHL2、FOSB昶WNT2B。

本研究對成年中華大蟾蜍耳后腺組織進行單細胞轉錄組測序，旨在鑒定蟾蜍耳后腺所有細胞，并對耳后腺細胞進行細胞亞群分類，從差異基因著手探索細胞異質性，篩選免疫基因，為提高中華大蟾蜍耳后腺取酥后存活率提供可能，為豐富蟾酥資源提供參考。

1材料與方法

1.1試驗材料

本研究所用中華大蟾蜍采自山東康宇中藥材有限公司中華大蟾蜍養殖基地，經山東中醫藥大學張永清教授鑒定為蟾蜍科蟾蜍屬的兩棲動物中華大蟾蜍（Bufo bufo gargarizans Cantor）。于2021年8月對未取酥和取酥后0、1、2、4、8、12、24、48、72、120、168 h的中華大蟾蜍進行耳后腺取樣，液氮速凍后-80℃冰箱保存備用。

1.2 scRNA-seq文庫構建

委托GeneDenovo生物技術有限公司（中國廣州），利用Illumina測序平臺的PE150測序模式進行高通量測序。

1.3測序數據分析

利用10×Genomics官方分析軟件Cell Ranger對原始數據進行過濾、比對、定量，鑒定回收細胞，去除低質量的reads后對每個樣本的reads數和測序質量進行初步統計，測序結果上傳GSA數據庫（GSA編號：GRA012164）。

1.4細胞亞群分類

使用Seurat軟件進行細胞聚類分析，保留細胞基因數量在330～2300范圍內、有效分子條形碼（unique molecular identifier，UMls）總數小于5600且單細胞中線粒體基因表達量比例小于10.0%的細胞。利用Harmony算法對數據進行合并及批次效應矯正，對合并后的數據進行主成分分析（principal component analysis，PCA）降維，采用soft k-means clustering算法對降維后的數據進行聚類。基于細胞亞群分類的結果，進一步利用t-分布領域嵌入算法（t-distributed stochastic neighbor embedding，tSNE）對細胞亞群分類結果進行可視化分析。

1.5細胞亞群上調表達基因分析

采用Seurat的秩和檢驗分別對不同類細胞群進行基因差異表達分析，以目標亞群或對照亞群中基因在25%以上的細胞中有表達、P≤0.01、基因表達倍數lOg2FC≥0.36為指標，篩選上調表達基因。然后利用Gene Ontology（GO）和KEGG數據庫對各聚類中基因的功能富集情況進行分析，包括GO功能分類注釋、GO功能顯著性富集分析及在差異基因中顯著性富集的通路。根據GO和KEGG分析結果，篩選與蟾蜍耳后腺取酥后免疫修復相關的基因。

1.6實時熒光定量PCR分析

使用TB Green Premix Ex TaqⅡ和TB Green Fast Primx EX TaqⅡ試劑盒，對成年中華大蟾蜍未取酥和取酥后不同時間段的耳后腺RNA樣本進行反轉錄后進行qRT-PCR分析。選擇在蟾蜍生物體中穩定表達的β-actin作為qRT-PCR驗證的內參基因，引物采用Primer3軟件（ht-tps：//bioinfo.ut.ee/primer3-0.4.01）設計，由Bio-Sune（中國上海）合成。引物序列見表1。

2結果與分析

2.1測序數據分析

利用10×Genomics官方分析軟件Cell Ranger對原始數據進行初步分析，共捕獲8414個細胞用于文庫構建和對端測序，結果如表2所示。共獲得389788801個Reads，其中有效條形碼占92.3%，有效分子條形碼（unique molecular identifi-er，UMls）占99.9%，81.2%的Reads可定位于中華大蟾蜍基因組。濾除細胞基因數小于330和超過2300、UMls超過5600和線粒體基因在單細胞的表達超過10.0%的細胞后，7687個高質量細胞被保留并進行后續分析。

2.2細胞亞群分類

利用Seurat軟件的典型相關分析函數對7687個高質量細胞進行聚類分析，通過tSNE非線性聚類的方法得到12個聚類（圖1）。每個聚類中的細胞數從34～4053個不等，占細胞總數的0.46%～52.37%（圖2）。對獲得的12個聚類進行細胞亞群識別，結果（圖3）顯示，Cluster 0中有紅細胞相關基因RHD、Alas2的表達，初步推測為紅細胞：Clus-ter 1中COLIA1、SMOC2的表達水平高于其他亞群，推測為來源于結締組織的成纖維細胞：Cluster 3中類固醇通路相關基因，如HSD1784、AKRID1、cyp2k4、SULT3A1和Cyp2c29存在高表達，推測可能是類固醇分泌細胞：Cluster 4和Cluster 7中存在高表達的內皮細胞marker基因SELE、VWF和CAV1，推測來源于導管的上皮細胞；Cluster 10中T細胞、B細胞抗原受體信號轉導和B細胞發育相關基因PTPRC、BLNK顯著表達，推測為來源于結締組織的漿細胞。Cluster 6中也存在基因PTPRC、BLNK的表達，但不如Cluster 10顯著，因此，選擇通過Clus-ter 10篩選蟾蜍耳后腺免疫修復基因。

2.3細胞亞群上調表達基因分析

為了解各細胞亞群的分子表達特征，采用Seurat的秩和檢驗對12個亞群中的上調基因進行篩選，并根據GO和KEGG分析其潛在功能。在上調基因的最小細胞百分比大于25%、平均表達量的10g2FC≥0.36的條件下，確定每個聚類中上調表達基因數量（圖4），然后選擇每聚類中豐度前5的基因繪制熱圖（圖5）。結果顯示，Clus-ter 10中的上調基因數為217個，涉及POU2AF1、IRF4、Cd79b和IGHC4等高表達基因。其中，POU2A F1主要在B細胞中表達，調控B細胞的發育和穩態，還是肺氣腫的關鍵調節因子，通過調節B細胞穩態參與慢性阻塞性肺疾病的進展；IRF4是相對分子質量為52 kDa的一種轉錄因子，其表達受限于免疫細胞，參與B細胞、T細胞、漿細胞、巨噬細胞等多種細胞的分化。

2.4 Cluster 10的GO和KEGG分析

對Cluster 10的上調基因進行GO功能注釋與分析，圖6展示GO分析的前20個GO terms，均為描述基因參與的生物過程（biological process），主要涉及淋巴細胞活化（lymphocyte ac-tivation）、淋巴細胞分化（lymphocyte differentia-tion）、免疫系統過程（lmmune system process）、白細胞活化（leukocyte activation）和正向調節生物過程（positive regulation of biological process）等，這些生物過程均與免疫系統有關。

對Cluster 10的KEGG富集分析顯示，前20個途徑包括了代謝（ metabolism）、環境信息處理（environmental information processing）、有機系統（organismal systems）及人類疾病（human diseases）（圖7A），其中有機系統及人類疾病為主要途徑；KEGG氣泡圖（圖7B）前5個富集途徑為B細胞受體信號通路（B cell receptor signaling pathway）、Fc epsilon RI信號通路（Fc epsilon RI signaling pathway）、趨化因子信號通路（chemokine signaling pathway）、原發性免疫缺陷（primary immunodefi-ciency）和Fc γ R介導的吞噬作用（Fc gamma R-mediated phagocytosis）。吞噬作用是最重要的免疫保護機制之一，Fc γ R可以介導刺激信號進入細胞，引起下游信號的改變，并導致肌動蛋白聚合和吞噬。而B細胞可作為免疫效應細胞直接參與體液免疫應答。

上述分析表明Cluster 10可能在機體的免疫應答和免疫功能中發揮重要作用，這與細胞亞群的鑒定結果相一致。

2.5 Cluster 10主要差異基因qRT-PCR分析

對Cluster 10中的差異基因，采用2^-ΔΔCt法計算相對表達量，并采用t-檢驗分析差異顯著性，篩選出MYOIF和CYBB兩個基因，對其在采酥后不同時間段的中華大蟾蜍耳后腺組織中的表達進行qRT-PCR分析。結果（圖8）顯示，兩基因相對表達量在取酥后12 h內并未有顯著變化，在24～48 h之間表達水平逐步升高，之后緩慢下降，推測兩個基因與免疫修復具有正向調節的關系。

3討論與結論

當動物機體受到損傷或感染時，免疫系統會自發抵抗病原微生物的入侵，并調節機體穩定。蟾蜍皮膚濕潤，皮膚分泌物是其抵抗微生物侵染、保護其免受環境和致病性損傷的第一道防線。而當蟾蜍取酥后，耳后腺受到損傷，皮膚分泌功能受到影響，使其防御能力大大降低，極易受到各種病菌的感染導致死亡。蟾蜍耳后腺相對于其他組織在相應甾族化合物的合成代謝上具有明顯優勢，已在其分泌物中鑒定出多種抗菌肽（Buforins，Peroniins，Brevinins等）。本研究使用10×Genomics系統對蟾蜍耳后腺進行scRNA-seq，鑒定蟾蜍耳后腺中的所有細胞亞群類型，研究細胞異質性，確定免疫細胞亞群，并篩選蟾蜍耳后腺免疫修復相關基因。

在scRNA-seq結果中，Cluster 10中GO富集分析的前20個GO terms均為基因參與的生物過程，其中MYOIF和CYBB參與了大部分生物過程，包括免疫系統過程、白細胞活化、正向調節生物過程、免疫應答、細胞活化、免疫系統過程的調節等。在KEGG分析通過Pathway顯著性富集確定差異基因參與的最主要生化代謝途徑和信號轉導途徑，其中MYOIF參與人類疾病中致病性大腸桿菌感染通路，CYBB參與有機系統中的白細胞經內皮細胞遷移通路。

通過qRT-PCR技術對MYOIF和CYBB兩個基因在蟾蜍耳后腺采酥后不同時間的相對表達量進行驗證，發現在采酥12h后其表達量隨時間的延長逐漸增加，在48h時表達量最高，推測這兩個基因對蟾蜍機體具有正向的免疫調控作用。MYOIF是一種非常規的肌球蛋白，通過控制肌動蛋白動力學在促進巨噬細胞介導的吞噬作用中發揮作用。MYOIF通過將AP2A1/αTAT1乙酰轉移酶復合物募集到α-微管蛋白上，使α-微管蛋白乙酰化。乙酰化α-微管蛋白在關鍵信號分子Syk和CARD9從細胞膜向細胞質的易位中起關鍵作用，這是抗真菌信號轉導的重要步驟。通過降低MYOIF表達，可以減少由IgE交聯和Mas相關的G蛋白偶聯受體X2（MRGPRX2）刺激誘導的人肥大細胞脫顆粒，且MYOIF對于AKT激活和DRP1磷酸化至關重要，是肥大細胞脫顆粒的重要調節因子。MYOIF取代VAV1的C端SH3結構形成的融合基因VAV1-myolf的表達可誘導CD4+T細胞淋巴瘤，CD4+、CD69+細胞增加，誘導共刺激受體（ICOS）上調，IL-10分泌增加。CYBB編碼黃細胞色素b558（也稱為gp91phox或NOX2）的β-鏈，而黃細胞色素b558是吞噬細胞（如粒細胞、單核細胞和巨噬細胞）中煙酰胺腺嘌呤二核苷酸磷酸（NADPH）氧化酶復合物（phox）的必需元素。CYBB在多種免疫細胞中表達，對IFN-γ、LPS和TNF-α等炎癥細胞因子產生反應，是吞噬細胞殺微生物氧化酶系統的主要組成部分。通過CRISPR/cas9介導的小鼠CYBB敲除，使轉基因小鼠產生殺微生物氧化酶缺陷，對自發微生物感染的敏感性增加，包括葡萄球菌、克雷伯氏菌以及真菌白色念珠菌。

本研究基于scRNA-seq技術，對中華大蟾蜍耳后腺取酥后免疫修復基因進行篩選。通過細胞亞群分析，確定Cluster 10為免疫細胞亞群。通過GO和KEGG對Cluster 10上調表達基因的分析，鑒定出多個參與免疫作用的生物過程與通路，并對篩選出的候選基因進行qRT-PCR驗證。根據qRT-PCR分析結果推測MYOIF和CYBB基因為蟾蜍耳后腺免疫修復正向調控基因。本研究可為后續通過生物技術手段對蟾蜍養殖進行人工干預提供數據支持，為挖掘蟾蜍耳后腺免疫修復基因，研究蟾蜍耳后腺取酥后的修復機制，提高蟾蜍采酥后成活率提供理論依據，同時為其他種蟾蜍耳后腺單細胞組學研究提供參考信息。