Identification of circular RNAs and their alterations involved in developing maleXenopus laevis chronically exposed to atrazine

Identification of circular RNAs and their alterations involved in developing maleXenopus laevis chronically exposed to atrazine 

    期刊：Chemosphere   影响因子：4.208分

    发表单位：山东省职业卫生与职业病防治研究院

1 摘要

    阿特拉津（Atrazine，简称AZ）是一种外源激素，可影响两栖动物正常生长发育。作者以往研究发现，AZ剂量100ug/L会引起发育期雄性蟾蜍发育异常，引起与性腺相关基因表达改变。最近，有人发现circRNA在多种发育异常的组织中参与作用。但是否有circRNAs参与AZ处理过的蟾蜍并不了解。本研究，通过对对照组（n=3）和AZ处理过的蟾蜍（n=3）的睾丸组织进行测序，鉴定到68575个circRNAs。AZ处理过组织中405个circRNAs表达有差异，其中44个上调，361个下调。通过qPCR对2个上调，6个下调的circRNAs进行验证，发现结果与转录组测序结果一致。282差异表达circRNAs发挥miRNA海绵吸附作用。KEGG通路分析发现，miRNAs靶基因主要涉及到19个pathway，其中Wnt通路， 孕激素介导的卵细胞成熟通路（progesterone-mediated oocyte maturation）可能是很重要的通路，本研究第一次提供证据显示，AZ能够影响雄性蟾蜍的正常发育，并能够调控睾丸中相关的基因表达。

2 研究方法

2.1 样品准备

    NF分期47（卵孵化后13d）, 来源于相同亲本，混合性别的蝌蚪，分成不同组，每个组（n=160），分成8个小组。AZ用二甲亚砜（DMSO, 0.01%, Sigma, St. Louis., MO, USA）溶解，处理组暴露于AZ（pu rity of 97%, Sigma, St. Louis., MO, USA）180天，对照组蝌蚪，用0.01% DMSO处理。

按照（Sai et al., 2016）方法，收集发育中雄性蟾蜍的睾丸，三个对照，三个处理，使用Trizol（Invitrogen, Carlsbad, CA, USA）提取TotalRNA，采用NanoDrop® ND-2000 光谱仪评估RNA的浓度和纯度，变性凝胶电泳检测RNA的完整性。

2.2 测序文库构建及circRNA测序分析

    使用去除核糖体RNA的试剂盒TruSeq Stranded Total RNA Library Prep Kit（Illumina, SanDiego, CA, USA） 构建文库。BioAnalyzer 2100 system （Agilent Technologies, Santa Clara, CA, USA）评估文库的质量。文库解链成单链DNA分子，测序仪流动槽捕获测，PCR反应成簇，Illumina HiSeq 4000 进行150个循环反应生成PE150序列。下机成对的序列，采用Q30控制质量，使用cutadapt软件去除3’接头及低质量reads。Clean reads将用于circRNAs分析，使用Bowtie2软件将序列比对到参考基因组/转录组（NCBI xl_ref_Xenopus_laevis_v2），使用find_circ软件鉴定circRNAs。成环剪切点reads，count值使用整体mapped reads进行校正，并进行log2转换。t-test进行两组差异表达分析，fold changes ≥2.0，P-values ≤0.05的circRNAs作为差异表达的circRNAs.

2.3 靶基因预测及生信分析

    使用Targetscan和Miranda 软件对circRNAs 结合miRNA位点，及miRNA靶基因进行预测，并对差异表达的circRNAs相关的靶基因进行KEGG富集分析。

2.4 circRNAs验证

    使用qPCR（Q-RT-PCR, ViiA 7 Real-time PCR System, Applied Biosystems, Carlsbad, CA, USA）对挑选的8个circRNA，进行验证，totalRNA 用 RNase R酶，消化线性RNA，引物采用divergent primer设计如Fig1，circRNAs的选择参考其潜在的功能和倍数两个因素。qPCR时，用Invitrogen Superscript cDNA Synthesis kit（Invitrogen Corp., Carlsbad, CA, USA）将TotalRNA转化成cDNA。RNase R酶处理时，TotalRNA（2ug）一部分用RNase R（20 U/μl, Epicentre）处理，一部分不用RNase R处理，37 °C 20min，使用2-△△CT法，计算表达量。2-△△CT>1时上调，2-△△CT<1下调，基因GPI用于内参。qPCR采用三重复结果以 mean±standard deviation形式呈现。使用SPSS Statistics 18.0 software（IBM Corp., New York, NY, USA）进行统计分析，使用student t-test at进行两组比较，p<0.05认为有差异。

3 研究结果

3.1 蟾蜍circRNAs鉴定

    对照组三个样品Clean reads分别得到85457078、85154262、70110988条，匹配到参考基因组的reads数分别是45908497、44827626 、37089055。经过AZ处理的处理组得到的Clean reads数目是89071634、95827312、84718576，匹配到参考基因组的reads数目是46499156、52443972、41544671。一共鉴定到68575 个circRNAs。

3.2 差异表达分析

处理组与对照组相比，共有405个circRNAs表达有差异，44个上调，361个下调，火山图如Fig.2。热图展示差异表达的circRNAs如Fig.3，行代表单个circRNAs，列代表样品名字，上面的树图代表样品聚类，左边的部分是AZ处理组，右边是对照组。本研究原始数据存储在GEO数据库中，编号为GSE103526。

3.3 circRNAs的miRNA结合位点预测，并构建circRNA-miRNA网络

    对差异表达的circRNAs进行miRNA预测，发现282个circRNAs中含有miRNA结合位点，结合位点数量从1-21不等。对这282个circRNAs选取前5个结合位点最多的miRNA, 如Table 4。然后自由选择了5个circRNAs，分别是NC_030735.1:109227383-109290634-, NC_030724.1:105611901-105637401+, NC_030734.1:106398702-106449286+, NW_016694872.1:550852-632532-, NC_030735.1:34874500-34913276-。构建circRNA-miRNA网络。如 fig4。

3.4 circRNA相关靶基因KEGG富集分析

    基于circRNAs相关的靶基因进行KEGG富集分析，发现19个信号通路，Table3, 其中，enrichment score排序发现，Wnt信号通路排在最前，并包含21个miRNAs靶基因，Wnt通路图如Fig.5。前十个通路展示如Fig.6。

3.5 circRNA验证

    对8个circRNAs进行qPCR发现，2个circRNAs是上调的（NC_030738.1:39776310-39778799- 和NC_030738.1:21126033-21131815+），6个circRNAs是下调的（NC_030734.1:102788414-102788574-, NW_016694872.1:550852-632532-,

NC_030737.1:31865134-31952348-, NC_030735.1:109227383-109290634-, NC_030734.1:106398702-106449286 + and NC_030732.1:43721272-43728946-;）qPCR验证的结果与circRNA测序的结果趋势一致。RNase R酶消化实验显示，这些RNA确是环状的。如图Fig7.

4 讨论

    AZ是一个常用的除草剂，对生态系统有一定的危害。暴露在AZ的环境，可能是两栖动物种群数量下降的一个重要原因，AZ能够影响他们的生长，代谢，免疫系统以及血细胞过程（Langerveld et al., 2009）。我们之前的研究显示AZ的毒性与蟾蜍基因水平的改变相关。AZ对蟾蜍的生存产生影响，而且睾丸的体重和性腺指数均显著下降。另外，AZ的不同剂量（0.01, 1, 10, 100μg/L）均会引起睾丸的退化，尤其是发育期的小蟾蜍暴露在100μg/L时。结果表明AZ通过干扰相关基因的表达量改变，影响雄性发育期蟾蜍的生殖能力和免疫系统。此结果也促使我们去研究其背后更多的潜在分子机制。本研究，我们用qPCR验证了RNase R酶处理下，circRNA的表达，结果显示我们的方法有效的鉴定到了蟾蜍睾丸中circRNAs的种类。circRNAs是一类近些年发现的一类非编码RNA，并受到广泛关注。circRNAs表达丰富，稳定，并且进化保守(笔者对circRNAs基因组来源序列多物种保守性分析发现确实很保守，circRNA大多来源于外显子剪切形成，外显子是基因编码区，不同物种间保守程度很高)，在胞内分子调控网络和表观遗传中通过吸附miRNAs起重要的调控作用。Memczak et al. 发现，circRNAs在不同的组织和发育时期的表达不同。自从基于对基因组测序的方式发现circRNAs以来, 大量的circRNAs被报道出来，包括人，古细菌，线虫，小鼠等（Salzman et al., 2012; Danan et al., 2012; Memczak et al., 2013）。但，在两栖动物中还没有通过RNAseq技术报道过circRNAs，本研究，通过RNAseq，我们发现了68575个circRNAs。之前没有关于AZ处理蟾蜍睾丸组织中circRNAs的数据。本研究第一次提供证据，表明circRNAs在受到AZ处理的蟾蜍睾丸组织中有差异表达。之前有报道称circRNAs的表达量变化与疾病的发生有关（Lin et al., 2016），因而，我们可以推断，circRNAs的改变与蟾蜍睾丸退化有关。

    Hansen et al. 报道circRNAs的吸附作用是一个普通的现象。circRNAs通过与miRNAs的结合，调控miRNAs的靶基因。我们发现405个差异表达的circRNAs中有282个circRNAs有miRNAs结合位点，数量1-21个不等。蟾蜍中circRNAs作为miRNAs海绵吸附作用与之前在人、鸡、小麦等物种中报道的类似（Memczak et al., 2013; Zhang et al., 2016; Wang et al., 2016）。我们之前的研究发现基因表达量改变在AZ处理的蟾蜍睾丸退化中发挥重要作用。因此，我们推断这些改变的circRNAs通过miRNA海绵吸附作用，阻止miRNAs与其调控的靶基因结合，进而在发育异常组织中发挥作用。

    我们进一步分析了circRNAs的靶基因的功能，在KEGG通路中共注释到19个显著的通路，并发现在Wnt通路中有21个靶基因（mRNA）.Wnt通路会受到胞内多种信号分子的调控，进而引起多种细胞响应，包括基因转录，增殖，分化，迁移等（Kikuchi et al., 2009; Grumolato et al., 2010; Simons et al., 2012），Wnt通路的精确调控在正常发育和组织维持中是极其重要的（Schneider et al., Rudloff and Kemler, 2012）。因而，我们推测，AZ处理的蟾蜍睾丸退化中，与Wnt信号通路相关的差异circRNAs起到了很重要的作用。以往研究，在差异表达的基因中我们发现了6个信号通路。孕酮介导的卵母细胞成熟信号通路（Progesterone-mediated oocyte maturation pathway）就是其中一个。有趣的是，此通路在本研究中也被鉴定到。因而我们推断，该通路也与AZ处理蟾蜍睾丸异常发育有关。这些数据同样给出了一个启示，差异表达的circRNAs可能可以作为新的雄性生殖毒理的分子标志物。而且，本研究为进一步研究AZ处理蟾蜍睾丸退化中circRNA调控分子机制，提供了一个坚实的基础。并且，本研究有助于我们理解circRNAs的一般功能，以及在环境如AZ响应中的作用。

5 结论

    本研究，在蟾蜍睾丸组织中鉴定到68575个circRNAs。 circRNAs的表达谱显示circRNA的表达量改变与睾丸生殖毒理密切相关。405个circRNAs差异表达，其中44个上调，361个下调。282个差异表达的circRNAs发挥了miRNA海绵吸附作用。KEGG通路显示，这些circRNAs主要涉及到19个通路，其中Wnt通路，Progesterone-mediated oocyte maturation通路是两个非常重要 的通路。我们的研究有助于理解，AZ处理的雄性蟾蜍生殖毒理中潜在的的分子机制。

本研究思路简要整理：

    研究目的：雄性蟾蜍睾丸退化中circRNAs是否参与该过程

    实验方法：

    材料：AZ处理的发育期雄性蟾蜍睾丸及未经AZ处理的对照，测序样品为3对3设计。

    RNAseq测序策略：做circRNA测序通常有两种建库方式，一种是去核糖体rRNA，RNase R去线性RNA，链特异性建库的方式；另一种是去核糖体rRNA, 不去线性RNA，而对所有长链RNA，进行链特异性建库的方式。本研究采用第二种方式，没有去处线性RNA。采用HiSeq 4000上机，PE150测序模式。

    circRNA差异分析：先使用find_circ对clean reads进行处理，获得circRNA count值，再对count值，使用total mapped reads进行校正，并做log2转换。然后使用t-test方法进行差异分析。fold changes ≥2.0，P-values ≤0.05的circRNAs作为差异表达的 circRNAs.

    海绵吸附作用miRNA预测：使用Targetscan 和Miranda 软件预测circRNA结合的miRNA（没有miRNA的表达量数据，此miRNA会用该物种所有miRNA进行扫描），由于一个circRNA可能会与几十，几百个miRNA结合，此处选择前5个结合位点最多的miRNA，作为circRNA的结合miRNA（此处确定miRNA的目的是为了进行构建，circ-miRNA网络及miRNA靶基因预测）。

    靶基因富集分析：circRNA-associated target genes富集，通常有两类，一是circRNA的来源基因，另一类是circRNA 结合的miRNA，miRNA的靶基因。作者对circRNA潜在结合的miRNA进行了预测，并做了miRNA 的靶基因预测。KEGG分析，基于miRNA的靶基因。

    circRNA qPCR验证：挑选8个进行qPCR验证，采用RNase R处理的和没有处理的两种方式。405个有差异的，circRNA选择哪些呢，作者考虑circRNA的表达量变化差异（差异比较大），及circRNA的重要性（此处考虑circRNA 结合miRNAs进而调控的重要基因，选择重要的circRNA）。

    文章结果展示： 

    1. 对下机数据clean reads及mapped reads进行统计描述。

    2. 差异表达的circRNA进行了统计描述，并用火山图和热图展示。

    3. Circ-miRNA结合，选取了5个circRNA（可以结合miRNA的circRNA有282个，如果都用，这个网会非常的大），构建circ-miRNA网络。

    4. KEGG对circRNA 结合miRNA 的靶基因进行富集分析，富集到19个通路，附件表格列出19个通路的信息，正文用散点图展示top10个通路，作者重点关注 Wnt通路，并在文中放了Wnt通路图，标注颜色为miRNA的靶基因。

    5. circRNA的验证，qPCR. 用柱状图说明，circRNA鉴定的准确性。

微信扫描二维码关注公众号获取更多相关文献信息