海洋与湖沼  2023, Vol. 54 Issue (3): 786-798   PDF    
http://dx.doi.org/10.11693/hyhz20220900225
中国海洋湖沼学会主办。
0

文章信息

李博, 谢玉素, 张留所. 2023.
LI Bo, XIE Yu-Su, ZHANG Liu-Suo. 2023.
海洋线虫Litoditis marina低氧胁迫响应的转录组分析
TRANSCRIPTOME ANALYSIS OF THE RESPONSE OF MARINE NEMATODE LITODITIS MARINA TO HYPOXIA STRESS
海洋与湖沼, 54(3): 786-798
Oceanologia et Limnologia Sinica, 54(3): 786-798.
http://dx.doi.org/10.11693/hyhz20220900225

文章历史

收稿日期:2022-09-01
收修改稿日期:2022-10-24
海洋线虫Litoditis marina低氧胁迫响应的转录组分析
李博1,2,3,4, 谢玉素1,2,3, 张留所1,2,3     
1. 中国科学院海洋研究所 中国科学院实验海洋生物学重点实验室 山东青岛 266071;
2. 青岛海洋科学与技术试点国家实验室海洋生物学与生物技术功能实验室 山东青岛 266237;
3. 中国科学院海洋大科学研究中心 山东青岛 266071;
4. 中国科学院大学 北京 100049
摘要:全球气候变化导致海洋氧含量降低, 海水低氧胁迫对海洋无脊椎动物的生长、发育和繁殖造成了严重的威胁。以海洋线虫Litoditis marina为实验对象, 观察了其在不同氧浓度(21%、3%、1%和0.5%)条件下的生长发育速率, 并对不同氧浓度条件下的L1幼虫样品进行了转录组测序和分析。实验结果表明, 当环境氧浓度从21%下降至3%时, L. marina的发育成熟速度明显加快, 进一步下降至1%时, 产卵时间延长并和21%氧浓度接近, 但当氧浓度为0.5%时, L. marina的产卵时间显著延迟。比较转录组分析表明, 相比于21%氧浓度环境, 3%、1%和0.5%低氧条件下L. marina的糖酵解、糖异生、硫代谢和线粒体碳代谢等通路相关基因的表达显著上调; 而寿命调控通路、细胞色素P450代谢通路和ABC转运蛋白相关基因的表达显著下调。研究结果发现的海洋线虫应对氧浓度胁迫的基因表达变化模式, 为深入理解海洋无脊椎动物应答低氧胁迫分子机制提供了重要参考。
关键词海洋线虫    Litoditis marina    低氧胁迫    糖酵解    硫代谢    碳代谢通路    
TRANSCRIPTOME ANALYSIS OF THE RESPONSE OF MARINE NEMATODE LITODITIS MARINA TO HYPOXIA STRESS
LI Bo1,2,3,4, XIE Yu-Su1,2,3, ZHANG Liu-Suo1,2,3     
1. CAS Key Laboratory of Experimental Marine Biology, Institute of Oceanology, Chinese Academy of Sciences, Qingdao 266071, China;
2. Laboratory for Marine Biology and Biotechnology, Pilot National Laboratory for Marine Science and Technology (Qingdao), Qingdao 266237, China;
3. Center for Ocean Mega-Science, Chinese Academy of Sciences, Qingdao 266071, China;
4. University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: Global climate change has led to ocean deoxygenation, and the hypoxic stress poses a serious threat to the growth, development and reproduction of marine invertebrates. In this study, the marine nematode Litoditis marina was exploited to analyze the impact of hypoxia to marine invertebrates. The L. marina L1 larvae were sampled for transcriptome sequencing and analysis and the growth and development rate under different oxygen concentrations (21%, 3%, 1% and 0.5%) were observed. Results showed that when ambient oxygen concentration decreased from 21% to 3%, the egg-laying time of L. marina was significantly accelerated, and when it further decreased to 1%, the egg laying time of L. marina was similar to that in 21% oxygen. When the oxygen concentration dropped to 0.5%, the egg laying time was significantly attenuated. The KEGG enrichment analysis showed that genes in several pathways such as glycolysis, gluconeogenesis, sulfur metabolism, and mitochondrial carbon metabolism were significantly up-regulated in 3%, 1% and 0.5% oxygen conditions, compared to those in 21% oxygen. However, the expression levels of genes in certain pathways such as nematode lifespan regulation, cytochrome P450 metabolism, and ABC transporters were significantly down-regulated under hypoxia conditions. This study provided a foundation to further study the functions of key genes in response or adaptation to hypoxia stress, in the context of global climate change and ocean deoxygenation.
Key words: marine nematode    Litoditis marina    hypoxia stress    glycolysis    sulfur metabolism    carbon metabolism pathway    

人类活动导致温室效应加剧, 全球变暖日益严重。而升温则会造成海水溶氧量下降, 研究表明海洋低氧胁迫影响鲍(Haliotis rufescens)、海湾扇贝(Argopecten irradians)和硬壳蛤(Mercenaria mercenaria)等海洋无脊椎动物的生长发育和生存(Kim et al, 2013; Gobler et al, 2014)。转录组分析显示低氧应答启动珊瑚低氧诱导转录因子HIF-1及其靶基因表达, 热激蛋白HSP-70和HSP-90在珊瑚低氧耐受中发挥重要作用(Alderdice et al, 2021)。而严重缺氧显著提升双壳贝类Lembulus bicuspidatus无氧糖酵解、抗氧化和质量控制通路基因的表达(Amorim et al, 2021)。小头虫Capitella teleta通过低氧感知蛋白TRPA1识别低氧胁迫的严重程度进而决定是否逃离低氧环境(Ogino et al, 2019)。环节动物Hermodice carunculata在中度缺氧胁迫下仅上调HIF-1α亚基的表达水平, 在严重缺氧条件下HIF-1β亚基的表达水平显著上调, 而在间歇性缺氧环境下其血红素转运功能相关基因显著上调(Grimes et al, 2021)。

模式生物秀丽线虫Caenorhabditis elegans适应1%和0.5%低氧环境, 需要低氧诱导因子hif-1基因(Jiang et al, 2001; Padilla et al, 2002; Nystul et al, 2004; Powell-Coffman, 2010)。研究发现秀丽线虫在低氧胁迫环境下110个基因的表达发生了显著变化, 其中63个依赖于hif-1基因(Shen et al, 2005)。研究报道核激素受体NHR-49与低氧诱导因子HIF-1平行调控秀丽线虫低氧胁迫应答(Doering et al, 2022)。此外, 蛋白翻译和胰岛素信号通路也参与调控秀丽线虫的低氧胁迫响应(Anderson et al, 2009; Mabon et al, 2009; Itani et al, 2021; Mulroney et al, 2021; Hemphill et al, 2022)。ddx-52基因编码秀丽线虫RNA解旋酶参与核糖体RNA的加工, 该基因功能缺失型突变体具有低氧耐受表型, 研究表明核糖体和tRNA合成协同调控秀丽线虫低氧胁迫应答(Itani et al, 2021; Mulroney et al, 2021)。

海洋线虫Litoditis marina在世界各地的沿海和口岸广泛分布, 在海洋生态系统起着重要作用(Derycke et al, 2016; Xie et al, 2020; Zhao et al, 2021), 但L. marina响应低氧胁迫机制的相关研究尚无报道。本研究以L. marina为实验材料, 通过RNA-seq初步探讨了L. marina响应低氧胁迫的发育响应与转录组变化特征, 以期为进一步利用基因编辑等技术方法深入研究在全球气候变化背景下, 海洋无脊椎动物应答低氧胁迫的分子机制奠定基础。

1 材料与方法 1.1 海水NGM平板制备

海水NGM (Nematode Growth Medium)培养基制备参照Xie等(2021)。新制备的海水NGM大板(直径9 cm)于室温下放置48 h后接菌或者置于4 ℃冰箱冷藏保存。每个用于高密度培养维持海洋线虫的海水NGM大板接种15×浓缩的大肠杆菌OP50 950 mL并且用灭菌的涂布棒涂匀, 室温放置48 h, 之后将涂菌的海水NGM大板置于4 ℃冰箱冷藏备用; 用于氧浓度处理收样的大板则接种500 μL 1× OP50并用涂布棒涂匀, 室温放置2 h后直接使用或4 ℃冷藏保存。表型观察使用的海水NGM小板(直径3.5 cm), 制备后在室温放置两天后每板接种30 μL 1× OP50菌液, 室温放置2 h后置于4 ℃冷藏备用。所有在4 ℃冷藏保存的海水NGM平板使用前需要在室温放置30 min。

1.2 海洋线虫L.marina的大规模同步培养和收样

待NGM大板培养的L. marina大量产卵时, 用灭菌海水将培养基上的线虫与卵洗涤至50 mL离心管内, 自然沉降5 min, 转移底部沉降下来的成虫到新的大板上备用, 将含有虫卵的上悬液用玻璃管转移至另一个新的50 mL离心管中, 3 000 g离心2 min。

去上清, 保留底部沉淀和约1 mL的液体, 转移到新的15 mL离心管中, 加入灭菌水到12 mL, 震荡重悬待底部沉淀散匀后1 300 g离心1 min, 重复上述步骤3次。第三次洗涤完成后吸出上清液, 加入灭菌水直至3 mL, 再加入3 mL的bleach液后震荡90 s, 立刻加入灭菌水补到12 mL, 1 300 g离心1 min。吸出上清液, 加入灭菌水至12 mL, 震荡重悬后1 300 g离心1 min。重复上述步骤2次, 最后将灭菌水替换为灭菌海水洗涤离心一次。

吸出上清液, 保留沉淀的虫卵和少许海水, 将卵吹散后, 通过400目网筛过滤, 转移至含10 mL灭菌海水的中培养皿(直径6 cm)中孵育20~22 h。通过450目网筛过滤, 将L1幼虫转移到15× OP50的海水NGM大板培养。第5天海水大板上L. marina大量产卵, 可用于扩大化同步培养或者L1幼虫样品收集。当海洋线虫的培养规模到达收样所需规模时可以开始进行多批次的液体bleach以及收样。

1.3 海洋线虫L. marina的生长观察实验

经过液体bleach虫卵, 孵化20~22 h后, 通过450目网筛过滤, 将L1幼虫转移到接种有20 μL大肠杆菌1× OP50的海水NGM小板(直径3.5 cm), 每个小板转入70只同步化的新孵化L1幼虫。在20 ℃条件下, 实验设为21%、3%、1%和0.5%四个氧浓度, 每个氧浓度设置三个生物学重复, 分别记录每个小板上的产卵时间。采用氧气控制设备型号为BioSpherix C21, 生化培养箱型号为上海龙跃LBI-275-N。

1.4 RNA-seq测序及分析

RNA-seq转录组测序实验包括样品检测、文库构建、文库质控和上机测序, 具体步骤参照Xie等(2021)。转录组数据分析前对测序获得Reads进行质控, 去除含有接头和低质量的Reads。采用DESeq2进行样品组间的差异表达分析, 使用edgeR进行差异分析, 差异表达基因检测中将差异倍数(Fold Change)≥1.5以及错误发现率(False Discovery Rate, FDR) < 0.05作为筛选标准。随后进行KEGG和GO富集分析。

1.5 荧光定量PCR

本实验选取部分表达发生显著性差异的关键基因进行qPCR验证(pck-2pgk-1pygl-1mmcm-1ctl-2cysl-2sdhb-1mce-1), 参考基因为cdc-42表 1为相关引物序列, 实验步骤参照Cao等(2022)。荧光定量PCR基因表达倍数变化, 通过delta-delta Ct method计算各个实验组三个生物学样本平均值, 与相应基因在21%氧气浓度下的比值, 使用Graphpad Prism 9通过双尾t检验计算P值, 显著性标准为P < 0.05。荧光定量PCR结果与RNA-seq间的相关性通过Pearson correlation analysis计算。

表 1 qPCR验证基因的引物序列 Tab. 1 The primer sequences of verification genes in qPCR
基因名 正向引物 反向引物
EVM0013809_cdc-42 GCCACTAAGGCTATGTACGACCT TTCAACGTGAACCTTCTTTTCAA
EVM0015867_pck-2 GAAAGAGGCATGCTCTCTCCCCT ACACCGTCTGGTACTTGTCG
EVM0011074_pgk-1 GAGGATGCAACCAGCAGTGT AACGGTCTTTGCACGGAGCA
EVM0005891_pygl-1 CTCTCCATTCCGATCTTCTCA GATTGGAGAGAAGCAGCCAG
EVM0005533_mmcm-1 GTCGAGTCTGGCATGACGAA CAATCGATACTTGTTCACTCC
EVM0011182_ctl-2 GTGCTCACAATACGCAACGC CCTGAGGCATGGAATACGCT
EVM0007291_cysl-2 CTATGAGTGTGGAACGCAGA TCCTTGGCTAGTTCTTGAGCTC
EVM0008590_sdhb-1 ATCTCGGTCCTGCTGTACTCAT GCAGTTCATGATGGTGTGACAT
EVM0002958_mce-1 CAACTCGAAGCTAGAGCTCCT GTGGATGTCTTTCACCTCAA
注: cdc-42为参考基因, 其余为验证基因
2 结果与分析 2.1 海洋线虫L. marina在不同氧气浓度条件下的产卵时间比较

将同步孵化的70只L. marina L1幼虫转移到接种有20 μL OP50菌液的海水NGM小板上, 每个氧浓度均包括三次重复实验, 每次重复实验对象为70只幼虫。如图 1所示, 在21%氧浓度条件下, L. marina产卵的平均时间为87.3 h (86、86和90 h); 3%氧浓度条件下产卵时间均值为79.7 h (78、80和81 h); 1%氧浓度条件下产卵平均时间为88 h (86、88和90 h); 0.5%氧浓度条件下产卵的平均时间为97.3 h (96、98和98 h)。结果显示与21%氧浓度相比, 在3%氧浓度下, 线虫发育显著加速, 而在0.5%氧浓度下, 发育速率显著延迟。

图 1 L. marina在不同氧浓度下的生长速率 Fig. 1 The growth rate of L. marina under different oxygen concentrations 注: 以21%为对照组, 采用Dunnett校正的单因素方差分析(ANOVA with Dunnett’s correction)检验实验组和对照组之间的显著性, 误差条代表对照样本和重复实验的标准误差(*P < 0.05)
2.2 海洋线虫L. marina在不同氧气浓度条件下的转录组分析

通过RNA-seq对海洋线虫L. marina不同氧浓度条件下的L1幼虫进行了转录组测序和分析, 各对比组中上调和下调基因的数量如表 2所示。图 2是差异表达基因的聚类热图。图 3图 4分别是各对比组中显著上调和显著下调的KEGG富集通路。

表 2 每个对比组合中差异基因的数目 Tab. 2 The number of differentially expressed genes in each comparison set
氧浓度比较组 总差异基因数 上调基因数 下调基因数
3% vs 21% 931 434 497
1% vs 21% 1 759 1 143 616
0.5% vs 21% 827 416 411

图 2 L. marina在不同氧浓度条件下差异基因的聚类热图 Fig. 2 Heatmap of differentially expressed genes of L. marina under different oxygen concentrations 注: 横坐标为样品名, 纵坐标为差异基因FPKM归一化后的数值, 红色代表上调, 蓝色表示下调

图 3 显著上调基因的KEGG富集通路 Fig. 3 KEGG enrichment of significantly up-regulated pathways

图 4 显著下调基因的KEGG富集通路 Fig. 4 KEGG enrichment of significantly down-regulated pathways
2.3 低氧条件下糖酵解/糖异生通路相关基因的表达显著上调

通过KEGG富集分析发现, 相较于21%氧浓度环境, 糖酵解/糖异生通路多个基因(gpi-1aldo-2gpd-3pgk-1pyc-1pck-2fbp-1)在3%、1%和0.5%三个低氧浓度条件的表达发生了显著性上调(图 5a)。其中糖酵解通路基因gpi-1aldo-2gpd-3pgk-1分别编码葡萄糖-6磷酸异构酶(glucose-6-phosphate isomerase)、二磷酸果糖醛缩酶(Fructose bisphosphate aldolase)、甘油醛-3-磷酸脱氢酶(glyceraldehyde 3-phosphate dehydrogenase)和磷酸甘油酸激酶(phosphoglycerate kinase), 而糖异生通路基因pyc-1pck-2fbp-1则分别编码丙酮酸羧化酶(pyruvate carboxylase)、磷酸烯醇式丙酮酸羧基激酶(phosphoenolypyruvate carboxykinase)和果糖1,6二磷酸酶(fructose-1,6-biphosphatase)。

图 5 低氧胁迫环境下L. marina上调通路相关基因的表达 Fig. 5 The expression of up-regulated pathway related genes of L. marina under hypoxia 注: a. 糖酵解/糖异生相关基因的表达水平; b. 硫代谢通路相关基因的表达水平; c. 线粒体碳代谢相关基因的表达水平。表达倍数变化表示处理组(3%、1%和0.5%)与对照组(21%)的基因表达量的比值(实验组FPKM值/对照组FPKM值)。采用独立样本t-test方法检验, 误差条代表3次重复实验的标准误差(*P < 0.05, **P < 0.01, ***P < 0.001)
2.4 低氧条件下硫代谢通路相关基因的表达显著上调

KEGG富集分析发现, 相较于21%氧浓度, 硫代谢通路多个基因在3%、1%、0.5%低氧浓度下发生了显著上调, 包括ethe-1mpst-1cysl-2cysl-3 (图 5b)。其中ethe-1编码过硫化物加双氧酶(ETHE1 persulfide dioxygenase), mpst-1编码3-巯基丙酮酸硫转移酶(3-mercaptopyruvate transferase), cysl-2cysl-3则是半胱氨酸合酶(cysteine synthase)编码基因。

2.5 低氧条件下线粒体碳代谢相关基因的表达显著上调

通过KEGG富集分析, 发现相较于21%氧浓度, 碳代谢通路多个线粒体功能相关基因(mce-1mmcm-1sdhb-1)在3%、1%和0.5%低氧浓度处理条件下的表达发生了显著上调(图 5c)。其中mce-1mmcm-1sdhb-1分别编码甲基丙二酰辅酶A差向异构酶(methylmalonyl-CoA epimerase)、甲基丙二酸单酰辅酶A变位酶(methylmalonyl-CoA mutase)和琥珀酸脱氢酶B亚基(succinate dehydrogenase complex subunit B)。

2.6 低氧条件下寿命调控通路相关基因的表达显著下调

相较于21%氧浓度, 寿命调控通路相关基因如hsp-70ctl-2zip-3在3%、1%和0.5%氧浓度下发生了显著下调, 而daf-2基因在3%和0.5%氧浓度下发生了显著下调(图 6a)。其中hsp-70ctl-2zip-3daf-2分别编码热休克蛋白、过氧化氢酶、bZIP转录因子和类胰岛素生长因子受体。

图 6 低氧胁迫环境下L. marina下调通路相关基因的表达 Fig. 6 The expression of down-regulated pathway related genes of L. marina under hypoxia 注: a. 线虫寿命通路相关基因的表达水平; b. 细胞色素P450通路相关基因的表达水平; c. ABC转运蛋白相关基因的表达水平。表达倍数变化表示处理组(3%、1%和0.5%)与对照组(21%)的基因表达量的比值(实验组FPKM值/对照组FPKM值)。采用独立样本t-test方法检验, 误差条代表3次重复实验的标准误差(*P < 0.05, **P < 0.01, ***P < 0.001)
2.7 低氧条件下细胞色素P450代谢通路与ABC转运蛋白相关基因的表达显著下调

KEGG富集分析表明, 相较于21%氧浓度, 细胞色素P450代谢通路的相关基因(ugt-19ugt-22ugt-23gst-5)和ABC转运蛋白相关基因(pgp-9pgp-14abt-4)在3%、1%和0.5%氧浓度下发生了显著性下调(图 6b, 6c)。ugt编码UDP-葡糖醛酸基转移酶(UDP-glucuronosyltransferase), gst-5基因编码谷胱甘肽S转移酶(glutathione S-transferase)。pgp基因编码P-糖蛋白(P-glycoprotein), abt-4编码ABC转运蛋白ABCA3(ATP binding cassette subfamily A member 3)。

2.8 过氧化物酶体和蜕皮相关基因的表达在0.5%氧浓度下发生显著下调

KEGG富集分析表明, 相较于21%氧浓度处理组, 过氧化物酶体通路的相关基因(acs-5acs-17ctl-2)和蜕皮相关基因(acn-1ptr-10ptr-14nas-28nas-31nas-36)在0.5%氧浓度下发生了显著下调(图 7a, 7b)。acs基因编码酰基辅酶A合成酶(acyl-CoA synthetase)。acn-1编码血管紧张素转换酶(angiotensin converting enzyme), ptr编码pathed受体蛋白, nas编码线虫虾红素类似物蛋白(nematode astacins-like protein)。

图 7 L. marina过氧化物酶体和蜕皮通路相关基因的表达 Fig. 7 The expression of genes related to peroxisome metabolism and molting of L. marina 注: a. 过氧化物酶体相关基因的表达水平; b. 蜕皮通路相关基因的表达水平。表达倍数变化表示处理组(3%、1%和0.5%)与对照组(21%)的基因表达量的比值(实验组FPKM值/对照组FPKM值)。采用独立样本t-test方法检验, 误差条代表3次重复实验的标准误差(*P < 0.05, **P < 0.01, ***P < 0.001)
2.9 荧光定量PCR验证

图 8a所示, 本研究选取了糖酵解代谢相关的pgk-1基因、糖异生代谢相关的pck-2基因、淀粉与蔗糖代谢相关的pygl-1基因、硫代谢相关的cysl-2基因、线虫寿命调控相关的ctl-2基因以及线粒体代谢相关的mce-1mmcm-1sdhb-1基因对转录组数据进行了qPCR验证, 简单线性回归分析结果表明qPCR检测得到的基因表达趋势和RNA-seq的表达趋势高度一致(图 8b)。

图 8 qPCR验证RNA-seq数据 Fig. 8 Validation of the RNA-seq results using qPCR 注: a. qPCR检测结果; b. qPCR和转录组表达数据间的关联分析
3 讨论 3.1 糖酵解和糖异生通路相关基因可能参与调控海洋线虫L. marina的低氧胁迫环境响应

研究表明人体细胞在1%和5%低氧环境下, 糖酵解通路相关基因的表达显著上调(Jain et al, 2020)。与之对应的是, 本研究发现相较于21%氧浓度, 海洋线虫L. marina糖酵解通路的多个基因(gpi-1aldo-2gpd-3pgk-1)在3%、1%和0.5%低氧环境下的表达显著上调(图 5a)。此外有研究者发现RNAi敲低gpd-3基因的表达会对秀丽线虫低氧条件下的运动活力造成负面影响, 但敲降gpi-1pgk-1基因则没有影响(Mendenhall et al, 2006)。在低氧条件下, 多数真核生物会改变代谢策略如线粒体呼吸作用供能减少, 糖酵解代谢上调以维持体内ATP水平的稳定, 这对于低氧胁迫下的细胞供能和生存至关重要(Kierans et al, 2021)。另外, 相较于21%氧浓度环境, 糖异生通路多个基因(pyc-1pck-2fbp-1)在低氧环境下(3%、1%和0.5%)的表达显著上调(图 5a)。与本研究相一致, Shen等(2005)报道发现在低氧胁迫环境下秀丽线虫pyc-1基因的表达显著上调, 同时发现pyc-1表达上调不依赖于低氧应答转录因子HIF-1。PCK-2催化草酰乙酸转换为磷酸烯醇式丙酮酸, 是糖异生代谢途径的限速步骤(Goncalves et al, 2020), 研究发现, pck-2过表达促进秀丽线虫运动能力, 而若RNAi敲低pck-2的表达水平, 则会损害秀丽线虫老年个体的运动能力, 并缩短平均寿命(Onken et al, 2020)。FBP-1催化果糖1,6二磷酸的磷酸酶水解, RNAi敲低fbp-1基因的表达减弱秀丽线虫老年成体的运动能力并显著缩短其平均寿命(Onken et al, 2020)。作者推测海洋线虫L. marina糖异生途径通过非碳水化合物生成的葡萄糖, 在低氧条件下很可能为显著增强的糖酵解途径提供原料且通过糖酵解供能应答低氧胁迫。

3.2 硫代谢途径相关基因可能参与了海洋线虫L. marina低氧胁迫响应

本研究发现相较于21%氧浓度, ethe-1mpst-1cysl-2cysl-3基因在3%、1%和0.5%低氧环境下的表达显著上调(图 5b)。ethe-1基因编码的硫加双氧酶在氧气参与下将源于H2S的多硫化合物氧化生成硫代硫酸盐, 参与H2S氧化通路(Budde et al, 2011)。mpst-1基因编码3-巯基丙酮酸硫转移酶, 秀丽线虫mpst-1功能缺失导致寿命缩减, 产卵量减少及H2S含量降低, 外源添加H2S可以营救寿命缩减表型(Qabazard et al, 2014)。cysl-2cysl-3编码硫化氢解酶, CYSL-2利用半胱氨酸和HCN催化合成H2S和β-氰基丙氨酸(Budde et al, 2011)。研究报道在低氧胁迫环境下秀丽线虫cysl-2基因的表达显著上调, 同时作者发现cysl-2基因的表达上调依赖于低氧应答转录因子HIF-1 (Shen et al, 2006)。与上述研究结果一致, 本研究发现与21%氧浓度相比, cysl-2基因在3%、1%和0.5%低氧浓度条件下的表达发生了显著上调(图 5b)。秀丽线虫cysl-2功能缺失突变体相比于野生型, 表现为产卵量减少, 内源性H2S含量降低, 寿命缩短(Qabazard et al, 2013)。H2S是一种具备重要生理学功能的生物体内源性小分子(Kabil et al, 2010)。内源性H2S作为一种抗氧化剂促进秀丽线虫寿命延长(Qabazard et al, 2014), 研究发现适量的H2S处理有助于激活秀丽线虫HIF-1因子(Budde et al, 2010)。作者推测3-巯基丙酮酸硫转移酶基因mpst-1和硫化氢解酶基因cysl-2cysl-3上调表达, 可能通过增加内源性H2S的合成以应答低氧胁迫。

3.3 线粒体碳代谢相关基因可能调控海洋线虫L. marina对低氧胁迫的响应

研究表明秀丽线虫mce-1突变体相比于野生型抗氧化胁迫能力显著增强(Kühnl et al, 2005)。mmcm-1编码蛋白催化R-甲基丙二酸单酰辅酶A转换为琥珀酰辅酶A (Bito et al, 2016)。琥珀酰辅酶A是三羧酸循环的重要组成环节, 在琥珀酰辅酶A合成酶的催化下, 通过高能酸酯键水解反应, 释放大量化学能从而引导GTP的合成。GTP可以和ATP互相转化, 在能量供给和细胞信号传导中都发挥了重要作用。并且琥珀酰辅酶A也是血红素合成的基本原料之一, 作为血红蛋白的组成成分, 血红素和血液氧气运输息息相关。sdhb-1基因编码琥珀酸脱氢酶B亚基, 是线粒体呼吸链中复合体Ⅱ的组成部分, 位于线粒体内膜, 参与TCA循环以及线粒体呼吸链的电子传递(Rustin et al, 2002)。秀丽线虫sdhb-1功能缺失突变体寿命缩短, 发育迟缓, 摄氧量降低, ATP合成减少(Saskői et al, 2020)。本研究发现相较于21%氧浓度, mce-1mmcm-1sdhb-1基因在3%、1%和0.5%三个低氧浓度条件下的表达显著上调(图 5c), 推测线粒体碳代谢相关基因可能参与了调控海洋线虫对低氧胁迫的应答。

3.4 线虫寿命调控相关基因参与了海洋线虫L. marina对低氧胁迫的响应

KEGG富集分析发现, 相较于21%氧浓度, 线虫寿命调控相关基因(ctl-2hsp-70daf-2zip-3)在3%、1%和0.5%三个低氧浓度条件下的表达发生了显著下调(图 6a)。ctl-2编码氧化氢酶, 其功能是将具有细胞毒性过氧化氢分解为无毒无害的成分, 主要作用部位是过氧化物酶体(Togo et al, 2000)。研究表明RNAi敲低ctl-2基因的表达导致秀丽线虫低氧响应因子HIF-1蛋白在体内积累(Xie et al, 2012), 据此推测CTL-2可能通过低氧胁迫应答因子HIF调控海洋线虫低氧环境应答。HSP70属于分子伴侣蛋白, 具有协助细胞内大分子组装和蛋白质折叠的功能, 其表达受多种环境胁迫的影响(Wieten et al, 2007)。研究发现在低氧胁迫环境下秀丽线虫hsp-70基因的表达显著上调, 同时发现hsp-70基因的表达上调不依赖于低氧应答转录因子HIF-1 (Shen et al, 2006)。与秀丽线虫低氧胁迫研究结果相反, 本研究发现海洋线虫hsp-70基因的表达在低氧胁迫环境下显著下调, 具体机制有待进一步研究。daf-2编码胰岛素受体(Kimura et al, 1997)。daf-2突变体增强了秀丽线虫对环境胁迫的抵抗能力, 并显著延长线虫寿命(Kenyon et al, 1993)。研究报道秀丽线虫daf-2突变体成体可以在20 ℃缺氧环境处理5 d后, 仍保持较高的存活率(Mendenhall et al, 2006)。与已经报道的秀丽线虫daf-2基因功能一致, 海洋线虫daf-2基因在低氧环境下显著下调, 表明胰岛素通路可能调控海洋线虫L. marina的低氧胁迫应答。zip-3基因编码转录因子, 秀丽线虫zip-3基因的下调表达激活线粒体未折叠蛋白反应(mitochondrial unfolded protein response UPRmt) (Deng et al, 2019)。UPRmt是线粒体在外界环境因子或者病原体胁迫下, 由于线粒体内变性以及错误折叠蛋白质积累, 所激发的由核DNA编码的线粒体应激反应, 目的是维持蛋白质稳态(Haynes et al, 2013)。因此作者推测在低氧胁迫环境下, 海洋线虫L. marina体内可能通过下调zip-3基因的表达维持蛋白质稳态。研究发现糖酵解抑制寿命和胁迫调控关键转录因子DAF-16并缩短秀丽线虫寿命(Onken et al, 2020), 本研究发现的寿命调控相关基因可能与糖酵解基因显著上调偶联。

3.5 细胞色素P450和ABC转运相关基因参与了海洋线虫L. marina对低氧胁迫的响应

细胞色素P450代谢通路的ugt-19ugt-22ugt-23gst-5基因在低氧浓度条件下发生了显著下调, 研究报道RNAi敲低ugt-19的表达促进秀丽线虫在低氧胁迫环境下的生存能力(Ladage et al, 2016)。RNAi敲低gst-5基因的表达导致秀丽线虫寿命缩短(Ayyadevara et al, 2007)。研究发现在低氧胁迫环境下秀丽线虫ugt-8基因的表达显著上调, 同时发现ugt-8基因的表达上调不依赖于低氧应答转录因子HIF-1 (Shen et al, 2006), 与此研究结果相反, 本研究发现海洋线虫在低氧环境下多个ugt基因的表达显著下调, ugt基因在低氧胁迫应答中的功能有待进一步研究。ABC转运蛋白相关的pgp-9pgp-14基因编码产物P-糖蛋白属于ATP依赖性的跨膜外排泵, 可以将多种大分子物质排出细胞(Raza et al, 2016)。abt-4基因编码的ATP依赖性转运蛋白参与脂质的跨膜运输(García-Montelongo et al, 2021)。作者推测在低氧胁迫环境下, 海洋线虫L. marina可能通过下调细胞色素P450代谢通路和ABC转运蛋白相关基因的表达调控特定大分子或脂类的生物合成和跨膜运输以应答低氧胁迫环境。

3.6 海洋线虫L. marina低氧胁迫响应与发育调控间的关联

本研究发现, 与21%氧浓度相比, 在3%氧浓度下, 线虫发育显著加速, 而在0.5%氧浓度下, 发育显著延迟。通过KEGG富集分析(p value=0.004 5), 相较于21%氧浓度, 过氧化物酶体通路相关基因(acs-5acs-17ctl-2)在0.5%氧浓度下发生了显著下调(图 7a)。ACS-5与ACS-17属于酰基辅酶A合成酶, 参与长链脂肪酸的代谢过程(Ruiz et al, 2019)。研究发现RNAi敲低acs-5基因抑制秀丽线虫由血清素处理引起的耗氧量增加(Srinivasan et al, 2008)。而RNAi敲低acs-17基因导致秀丽线虫产卵量下降, 胚胎致死率上升(Tischler et al, 2006)。与21%氧浓度相比, 在3%氧浓度下acs-5acs-17ctl-2基因表达没有显著变化(图 7a)。我们推测上述过氧化物酶体相关基因的表达下调可能是0.5%氧浓度下海洋线虫L. marina生长缓慢的主要原因之一。

通过GO富集分析(p value=0.000 2)发现, 相较于21%氧浓度, 蜕皮通路相关基因(acn-1ptr-10ptr-14nas-28nas-31nas-36)在0.5%氧浓度下发生了显著下调(图 7b)。研究表明acn-1基因的表达对于秀丽线虫幼虫蜕皮是必需的(Brooks et al, 2003), RNAi敲低acn-1基因的表达导致秀丽线虫幼虫发育缓慢(Kamath et al, 2003)。ptr-10ptr-14基因表达产物属于patched受体蛋白, RNAi敲低ptr-10导致幼虫蜕皮缺陷, 并抑制生长发育(Tischler et al, 2006)。ptr-14基因的RNAi敲低导致秀丽线虫生长过程中运动能力变差, 并抑制生长发育(Zugasti et al, 2005)。nas-28, nas-31nas-36基因编码线虫虾红蛋白酶(nematode astacins-like protease), 属于金属内肽酶(Möhrlen et al, 2003)。研究报道nas-36基因的RNAi敲低导致秀丽线虫蜕皮缺陷, 虫体后端难以蜕皮(Suzuki et al, 2004)。与21%氧浓度相比, 在3%氧浓度下acn-1ptr-14基因表达没有发生显著变化(图 7b), 虽然ptr-10nas-28nas-31nas-36等基因在3%氧浓度下表达显著下调, 但这些基因总体上在0.5%氧浓度下表达下调幅度更大(图 7b), 作者推测上述蜕皮相关基因的表达下调可能是0.5%氧浓度下海洋线虫L. marina生长缓慢的主要原因之一。

4 结论

本研究结果表明, 相较于21%氧浓度, 海洋线虫在低氧环境下糖酵解、糖异生、硫代谢和线粒体碳代谢等通路相关的多个基因的表达显著上调; 而线虫寿命调控、细胞色素P450和ABC转运蛋白等通路相关基因的表达在低氧胁迫条件下显著下调。海洋线虫L. marina中发现的低氧胁迫应答调控模式, 为理解海洋无脊椎动物适应低氧环境的分子机制提供了参考, 同时为筛选和验证生物体响应和适应低氧胁迫的关键基因奠定了重要基础。

参考文献
ALDERDICE R, SUGGETT D J, CÁRDENAS A, et al, 2021. Divergent expression of hypoxia response systems under deoxygenation in reef-forming corals aligns with bleaching susceptibility. Global Change Biology, 27(2): 312-326 DOI:10.1111/gcb.15436
AMORIM K, PIONTKIVSKA H, ZETTLER M L, et al, 2021. Transcriptional response of key metabolic and stress response genes of a nuculanid bivalve, Lembulus bicuspidatus from an oxygen minimum zone exposed to hypoxia-reoxygenation. Comparative Biochemistry and Physiology Part B: Biochemistry and Molecular Biology, 256: 110617 DOI:10.1016/j.cbpb.2021.110617
ANDERSON L L, MAO X R, SCOTT B A, et al, 2009. Survival from hypoxia in C. elegans by inactivation of aminoacyl-tRNA synthetases. Science, 323(5914): 630-633 DOI:10.1126/science.1166175
AYYADEVARA S, DANDAPAT A, SINGH S P, et al, 2007. Life span and stress resistance of Caenorhabditis elegans are differentially affected by glutathione transferases metabolizing 4-hydroxynon-2-enal. Mechanisms of Ageing and Development, 128(2): 196-205 DOI:10.1016/j.mad.2006.11.025
BITO T, WATANABE F, 2016. Biochemistry, function, and deficiency of vitamin B12 in Caenorhabditis elegans. Experimental Biology and Medicine, 241(15): 1663-1668 DOI:10.1177/1535370216662713
BROOKS D R, APPLEFORD P J, MURRAY L, et al, 2003. An essential role in molting and morphogenesis of Caenorhabditis elegans for ACN-1, a novel member of the angiotensin-converting enzyme family that lacks a metallopeptidase active site. Journal of Biological Chemistry, 278(52): 52340-52346 DOI:10.1074/jbc.M308858200
BUDDE M W, ROTH M B, 2010. Hydrogen sulfide increases hypoxia-inducible factor-1 activity independently of von Hippel-Lindau tumor suppressor-1 in C. elegans. Molecular Biology of the Cell, 21(1): 212-217 DOI:10.1091/mbc.e09-03-0199
BUDDE M W, ROTH M B, 2011. The response of Caenorhabditis elegans to hydrogen sulfide and hydrogen cyanide. Genetics, 189(2): 521-532 DOI:10.1534/genetics.111.129841
CAO X W, SUN P Q, ZHANG L S, 2022. Transcriptome analysis of the marine nematode Litoditis marina in a chemically defined food environment with stearic acid supplementation. Journal of Marine Science and Engineering, 10(3): 428 DOI:10.3390/jmse10030428
DENG P, UMA NARESH N, DU Y G, et al, 2019. Mitochondrial UPR repression during Pseudomonas aeruginosa infection requires the bZIP protein ZIP-3. Proceedings of the National Academy of Sciences of the United States of America, 116(13): 6146-6151 DOI:10.1073/pnas.1817259116
DERYCKE S, DE MEESTER N, RIGAUX A, et al, 2016. Coexisting cryptic species of the Litoditis marina complex (Nematoda) show differential resource use and have distinct microbiomes with high intraspecific variability. Molecular Ecology, 25(9): 2093-2110 DOI:10.1111/mec.13597
DOERING K R, CHENG X J, MILBURN L, et al, 2022. Nuclear hormone receptor NHR-49 acts in parallel with HIF-1 to promote hypoxia adaptation in Caenorhabditis elegans. eLife, 11: e67911 DOI:10.7554/eLife.67911
GARCÍA-MONTELONGO M, GONZÁLEZ-VILLARREAL S E, DEL RINCÓN-CASTRO M C, et al, 2021. Use of RNAi as a preliminary tool for screening putative receptors of nematicidal toxins from Bacillus thuringiensis. Archives of Microbiology, 203(4): 1649-1656 DOI:10.1007/s00203-020-02179-1
GOBLER C J, DEPASQUALE E L, GRIFFITH A W, et al, 2014. Hypoxia and acidification have additive and synergistic negative effects on the growth, survival, and metamorphosis of early life stage bivalves. PLoS One, 9(1): e83648 DOI:10.1371/journal.pone.0083648
GONCALVES J, WAN Y F, GUO X Y, et al, 2020. Succinate dehydrogenase-regulated phosphoenolpyruvate carboxykinase sustains copulation fitness in aging C. elegans males. iScience, 23(4): 100990 DOI:10.1016/j.isci.2020.100990
GRIMES C J, PETERSEN L H, SCHULZE A, 2021. Differential gene expression indicates modulated responses to chronic and intermittent hypoxia in corallivorous fireworms (Hermodice carunculata). Scientific Reports, 11(1): 11110 DOI:10.1038/s41598-021-90540-9
HAYNES C M, FIORESE C J, LIN Y F, 2013. Evaluating and responding to mitochondrial dysfunction: the mitochondrial unfolded-protein response and beyond. Trends in Cell Biology, 23(7): 311-318 DOI:10.1016/j.tcb.2013.02.002
HEMPHILL C, PYLARINOU-SINCLAIR E, ITANI O, et al, 2022. Daf-16 mediated repression of cytosolic ribosomal protein genes facilitates a hypoxia sensitive to hypoxia resistant transformation in long-lived germline mutants. PLoS Genetics, 18(5): e1009672 DOI:10.1371/journal.pgen.1009672
ITANI O A, ZHONG X F, TANG X T, et al, 2021. Coordinate regulation of ribosome and tRNA biogenesis controls hypoxic injury and translation. Current Biology, 31(1): 128-137 DOI:10.1016/j.cub.2020.10.001
JAIN I H, CALVO S E, MARKHARD A L, et al, 2020. Genetic screen for cell fitness in high or low oxygen highlights mitochondrial and lipid metabolism. Cell, 181(3): 716-727 DOI:10.1016/j.cell.2020.03.029
JIANG H Q, GUO R, POWELL-COFFMAN J A, 2001. The Caenorhabditis elegans hif-1 gene encodes a bHLH-PAS protein that is required for adaptation to hypoxia. Proceedings of the National Academy of Sciences of the United States of America, 98(14): 7916-7921 DOI:10.1073/pnas.141234698
KABIL O, BANERJEE R, 2010. Redox biochemistry of hydrogen sulfide. Journal of Biological Chemistry, 285(29): 21903-21907 DOI:10.1074/jbc.R110.128363
KAMATH R S, FRASER A G, DONG Y, et al, 2003. Systematic functional analysis of the Caenorhabditis elegans genome using RNAi. Nature, 421(6920): 231-237 DOI:10.1038/nature01278
KENYON C, CHANG J, GENSCH E, et al, 1993. A C. elegans mutant that lives twice as long as wild type. Nature, 366(6454): 461-464 DOI:10.1038/366461a0
KIERANS S J, TAYLOR C T, 2021. Regulation of glycolysis by the hypoxia-inducible factor (HIF): implications for cellular physiology. The Journal of Physiology, 599(1): 23-37 DOI:10.1113/JP280572
KIM T W, BARRY J P, MICHELI F, 2013. The effects of intermittent exposure to low-pH and low-oxygen conditions on survival and growth of juvenile red abalone. Biogeosciences, 10(11): 7255-7262 DOI:10.5194/bg-10-7255-2013
KIMURA K D, TISSENBAUM H A, LIU Y X, et al, 1997. daf-2, an insulin receptor-like gene that regulates longevity and diapause in Caenorhabditis elegans. Science, 277(5328): 942-946 DOI:10.1126/science.277.5328.942
KÜHNL J, BOBIK T, PROCTER J B, et al, 2005. Functional analysis of the methylmalonyl-CoA epimerase from Caenorhabditis elegans. The FEBS Journal, 272(6): 1465-1477 DOI:10.1111/j.1742-4658.2005.04579.x
LADAGE M L, KING S D, BURKS D J, et al, 2016. Glucose or altered ceramide biosynthesis mediate oxygen deprivation sensitivity through novel pathways revealed by transcriptome analysis in Caenorhabditis elegans. G3 Genes|Genomes|Genetics, 6(10): 3149-3160
MABON M E, SCOTT B A, CROWDER C M, 2009. Divergent mechanisms controlling hypoxic sensitivity and lifespan by the DAF-2/insulin/IGF-receptor pathway. PLoS One, 4(11): e7937 DOI:10.1371/journal.pone.0007937
MENDENHALL A R, LARUE B, PADILLA P A, 2006. Glyceraldehyde-3-phosphate dehydrogenase mediates anoxia response and survival in Caenorhabditis elegans. Genetics, 174(3): 1173-1187 DOI:10.1534/genetics.106.061390
MÖHRLEN F, HUTTER H, ZWILLING R, 2003. The astacin protein family in Caenorhabditis elegans. European Journal of Biochemistry, 270(24): 4909-4920 DOI:10.1046/j.1432-1033.2003.03891.x
MULRONEY T E, PÖYRY T, WILLIS A E, 2021. Hypoxia: uncharged tRNA to the rescue!. Current Biology, 31(1): R25-R27 DOI:10.1016/j.cub.2020.10.067
NYSTUL T G, ROTH M B, 2004. Carbon monoxide-induced suspended animation protects against hypoxic damage in Caenorhabditis elegans. Proceedings of the National Academy of Sciences of the United States of America, 101(24): 9133-9136 DOI:10.1073/pnas.0403312101
OGINO T, TOYOHARA H, 2019. Identification of possible hypoxia sensor for behavioral responses in a marine annelid, Capitella teleta. Biology Open, 8(3): bio037630
ONKEN B, KALINAVA N, DRISCOLL M, 2020. Gluconeogenesis and PEPCK are critical components of healthy aging and dietary restriction life extension. PLoS Genetics, 16(8): e1008982 DOI:10.1371/journal.pgen.1008982
PADILLA P A, NYSTUL T G, ZAGER R A, et al, 2002. Dephosphorylation of cell cycle-regulated proteins correlates with anoxia-induced suspended animation in Caenorhabditis elegans. Molecular Biology of the Cell, 13(5): 1473-1483 DOI:10.1091/mbc.01-12-0594
POWELL-COFFMAN J A, 2010. Hypoxia signaling and resistance in C. elegans. Trends in Endocrinology & Metabolism, 21(7): 435-440
QABAZARD B, AHMED S, LI L, et al, 2013. C. elegans aging is modulated by hydrogen sulfide and the sulfhydrylase/ cysteine synthase cysl-2. PLoS One, 8(11): e80135 DOI:10.1371/journal.pone.0080135
QABAZARD B, LI L, GRUBER J, et al, 2014. Hydrogen sulfide is an endogenous regulator of aging in Caenorhabditis elegans. Antioxidants & Redox Signaling, 20(16): 2621-2630
RAZA A, BAGNALL N H, JABBAR A, et al, 2016. Increased expression of ATP binding cassette transporter genes following exposure of Haemonchus contortus larvae to a high concentration of monepantel in vitro. Parasites & Vectors, 9(1): 522
RUIZ M, BODHICHARLA R, STÅHLMAN M, et al, 2019. Evolutionarily conserved long-chain Acyl-CoA synthetases regulate membrane composition and fluidity. eLife, 8: e47733 DOI:10.7554/eLife.47733
RUSTIN P, RÖTIG A, 2002. Inborn errors of complex Ⅱ-unusual human mitochondrial diseases. Biochimica et Biophysica Acta (BBA) - Bioenergetics, 1553(1/2): 117-122
SASKŐI É, HUJBER Z, NYÍRŐ G, et al, 2020. The SDHB Arg230His mutation causing familial paraganglioma alters glycolysis in a new Caenorhabditis elegans model. Disease Models & Mechanisms, 13(10): dmm044925
SHEN C, NETTLETON D, JIANG M, et al, 2005. Roles of the HIF-1 hypoxia-inducible factor during hypoxia response in Caenorhabditis elegans. Journal of Biological Chemistry, 280(21): 20580-20588 DOI:10.1074/jbc.M501894200
SHEN C, SHAO Z, Powell-COFFMAN J A, 2006. The Caenorhabditis elegans rhy-1 gene inhibits HIF-1 hypoxia-inducible factor activity in a negative feedback loop that does not include vhl-1. Genetics, 174(3): 1205-1214 DOI:10.1534/genetics.106.063594
SRINIVASAN S, SADEGH L, ELLE I C, et al, 2008. Serotonin regulates C. elegans fat and feeding through independent molecular mechanisms. Cell Metabolism, 7: 533-544 DOI:10.1016/j.cmet.2008.04.012
SUZUKI M, SAGOH N, IWASAKI H, et al, 2004. Metalloproteases with EGF, CUB, and thrombospondin-1 domains function in molting of Caenorhabditis elegans. Biological Chemistry, 385(6): 565-568 DOI:10.1515/BC.2004.069
TISCHLER J, LEHNER B, CHEN N S, et al, 2006. Combinatorial RNA interference in Caenorhabditis elegans reveals that redundancy between gene duplicates can be maintained for more than 80 million years of evolution. Genome Biology, 7(8): R69 DOI:10.1186/gb-2006-7-8-r69
TOGO S H, MAEBUCHI M, YOKOTA S, et al, 2000. Immunological detection of alkaline-diaminobenzidine-negativeperoxisomes of the nematode Caenorhabditis elegans: purification and unique pH optima of peroxisomal catalase. European Journal of Biochemistry, 267(5): 1307-1312 DOI:10.1046/j.1432-1327.2000.01091.x
WIETEN L, BROERE F, VAN DER ZEE R, et al, 2007. Cell stress induced HSP are targets of regulatory T cells: A role for HSP inducing compounds as anti-inflammatory immune-modulators?. FEBS Letters, 581(19): 3716-3722 DOI:10.1016/j.febslet.2007.04.082
XIE M, ROY R, 2012. Increased levels of hydrogen peroxide induce a HIF-1-dependent modification of lipid metabolism in AMPK compromised C. elegans dauer larvae. Cell Metabolism, 16(3): 322-335 DOI:10.1016/j.cmet.2012.07.016
XIE Y S, ZHANG P C, XUE B N, et al, 2020. Establishment of a marine nematode model for animal functional genomics, environmental adaptation and developmental evolution. bioRxiv DOI:10.1101/2020.03.06.980219
XIE Y S, ZHANG P C, ZHANG L S, 2021. Genome-wide transcriptional responses of marine nematode Litoditis marina to hyposaline and hypersaline stresses. Frontiers in Physiology, 12: 672099 DOI:10.3389/fphys.2021.672099
ZHAO L, GAO F, GAO S, et al, 2021. Biodiversity-based development and evolution: the emerging research systems in model and non-model organisms. Science China Life Sciences, 64(8): 1236-1280 DOI:10.1007/s11427-020-1915-y
ZUGASTI O, RAJAN J, KUWABARA P E, 2005. The function and expansion of the patched- and hedgehog-related homologs in C. elegans. Genome Research, 15(10): 1402-1410 DOI:10.1101/gr.3935405