干旱胁迫下桦木梨circRNAs的鉴定与表征

摘要

环状RNA（circRNA）在miRNA功能和转录控制中起重要作用。然而，关于梨中的circRNAs知之甚少。在这项研究中，我们使用深度测序鉴定了 circRNA，并分析了它们在干旱胁迫下的表达。我们总共鉴定了 899 个 circRNA，其中 33 个（23 个上调，10 个下调）显示出脱水反应。我们进行了 GO 和 KEGG 富集分析来预测差异表达的 circRNA 的功能。预计 309 个 circRNA 可充当 180 个 miRNA 的海绵。基于差异表达的circRNA与其miRNA结合位点之间的相关性分析，构建了circRNA-miRNA共表达 络。我们的研究将为发现与干旱胁迫相关的基因提供丰富的遗传资源，

介绍

根据最近的研究，真核基因组编码了大量的ncRNA（非编码RNA）[ 1 , 2 ]。ncRNA 没有或几乎没有蛋白质编码潜力，但在各种生物过程中发挥作用 [ 3 ]。环状 RNA (circRNA) 是一种 ncRNA 分子，它没有 5? 帽和聚 A 尾。circRNA是由共价键形成的环状结构，circRNA分别于1980年和1993年首先在酵母[ 4 ]和人类[ 5 ]中被发现。环状 RNA 是在剪接过程中通过各种机制产生的 [ 6 ]：大多数环状 RNA 是由反向剪接的外显子产生的，而其他的则来自内含子 [ 7]]。circRNA在细胞中比线性RNA分子更稳定，半衰期更长，可以抵抗RNAase R降解[ 8 ]。根据它们在基因组中的一般位置，有五种类型的 circRNA：外显子 circRNA、内含子 circRNA、有义重叠 circRNA、反义 circRNA 和基因间 circRNA [ 9 ]。

circRNA 在生命的所有领域中无处不在，可以实现多种生物学功能。它们可以作为竞争性内源 RNA 来结合微 RNA (miRNA) [ 10 , 11 ]、基因转录和表达调节剂 [ 12 ] 以及 RNA 结合蛋白（蛋白质海绵）[ 13 ]。例如，超过 70 个 miRNA 靶位点（选择性保守）存在于 ciRS-7 中，并且 ciRS-7 与人类 miR-7 依赖型的 Argonaute 蛋白密切相关 [ 10 ]。circRNA 在应激反应中也起着至关重要的作用。例如，与对照相比，163 个 circRNA 在暴露于冷应激的西红柿中显着差异表达 [ 14]。]。通过高通量测序平台和生物信息学方法的发展的辅助下，circRNAs丰盈最近已在古细菌鉴定[ 15 ]，人类，小鼠[ 16 ]，拟南芥[ 17，18 ]，和水稻[ 12 ]。

多项研究表明，压力会诱导多个circRNAs。在这项研究中，我们使用 RNA 测序 (RNA-seq) 分析了干旱胁迫下桦树叶梨叶中差异表达的 circRNA。我们的研究结果将促进梨育种计划的发展，并为鉴定其他果树中的 circRNA 奠定基础。

材料和方法

桦叶梨样品

桦叶梨（P . betulifolia Bunge ）幼苗在中国江苏省南京市江苏省农业科学院园艺研究所国家种质果园的苗床中生长。根据我们之前的 告 [ 20 ] ，将幼苗置于生长室中，以 24 小时为周期：25°C 光照 14 小时，20°C 黑暗 10 小时。将六叶期幼苗插入盛有蒸馏水的烧杯中 2 d，然后进行脱水处理。然后将幼苗转移到 1/2 Murashige 和 Skoog Basal 中(MS) 溶液含有 15% 聚乙二醇 (PEG) 以模拟干旱胁迫。在处理后 48 小时收集对照和处理幼苗的叶子，一式三份，用蒸馏水冲洗，在液氮中冷冻，并在 –80°C 下储存直至进一步使用。

RNA制备

使用 Takara Mini BEST 植物 RNA 提取试剂盒按照制造商的方案从六个样品中的每一个中提取总 RNA。使用甲醛琼脂糖凝胶电泳验证 RNA 质量，使用 NanoDrop ND-1000 分光光度计确定 RNA 数量。每个样品总共使用 5 μg RNA 用于实验。为了获得核糖体 RNA (rRNA) 去除的 RNA，使用 Epicenter Ribo-zero rRNA Removal Kit（Epicentre，美国）去除 rRNA。随后，按照制造商的建议，使用 NEBNext Ultra Directional RNA Library Prep Kit for Illumina（NEB，美国）从 rRNA 耗尽和 RNase R 消化的 RNA 生成测序文库。最后，使用 AMPure XP 系统纯化文库并使用 Agilent Bioanalyzer 2100 系统进行鉴定。

聚类和测序

根据制造商的说明，使用 HiSeq PE Cluster Kit v4 cBot (Illumina) 在 cBot Cluster Generation System 上对索引编码样本进行聚类。簇生成后，在 Illumina Hiseq 4000 平台上对文库进行测序，并生成 150 bp 双端读数。

质量控制

原始数据（fastq 格式的原始读取）首先使用自定义 Perl 脚本进行处理。从原始数据中去除包含适配器的读取、包含多聚 N 的读取和低质量读取后，获得干净数据（干净读取）。计算清洁数据的 Q20、Q30 和 GC 含量。所有下游分析均基于此步骤中生成的干净、高质量的数据。

映射到参考基因组

参考基因组和基因模型注释文件从梨基因组 站[ 21 ]（http://peargenome.njau.edu.cn/）下载。使用Bowtie v2.0.6软件构建参考基因组索引，使用TopHat v2.0.9软件将双端clean read与参考基因组比对。

circRNA鉴定

我们使用了 find_circ 版本 1.1 [ 11 ] 和 CIRI 版本 2.0.5 [ 22 ]] 来识别 circRNA。对于 find_circ，保留未映射的读数，并使用 Bowtie v2.0.6 提取这些读数的 5′ 和 3′ 末端的 20-mers 并与参考序列独立比对。锚定序列通过 find_circ 算法进行扩展，使得完整的读取对齐并且断点两侧是 GU/AG 剪接位点。接下来，将具有至少两个支持读数的反向拼接读数鉴定为 circRNA。CIRI 算法是另一种识别 circRNA 的工具，它扫描 SAM 文件两次并收集足够的信息来识别和表征 circRNA。CIRI 使用默认选项执行，识别的 circRNA 读数计数通过读数长度标准化，并在 CIRI 预测后确定读数映射数。最后，

量化

circRNA 表达水平基于每百万转录本 (TPM) 使用以下公式进行归一化 [ 23 ]：

归一化表达式= (映射读取) / (总计读取) ×1,000,000。

差异表达

使用DESeq2 1.6.3版软件[ 24 ]进行两组之间的差异表达（0h和48h）。使用 Benjamini 和 Hochberg 方法调整 P 值。默认情况下，差异表达的校正 p 值阈值设置为 0.05。

GO 和 KEGG 富集分析

使用GOseq软件（版本1.18.0）对差异circRNA的源基因进行基因本体论（GO，http://www.geneontology.org/）富集分析。KEGG [ 25 ] 是一个用于理解生物系统高级功能的数据库（http://www.genome.jp/kegg/）。KOBAS 络服务器 [ 26 ] 用于京都基因和基因组百科全书 (KEGG) 通路富集分析。

miRNA结合位点预测

CircRNA 可以通过与 miRNA 结合来抑制 miRNA 的功能。为了进一步研究circRNA的功能，我们使用PsRobot软件分析了鉴定的circRNA的miRNA结合位点。该软件是一个基于 络的易用工具，它使用 miRNA 和 circRNA 序列（包括连接位点）进行目标基因预测，它可以执行快速分析以使用批量输入数据中识别具有茎环形状前体的 miRNA。改进的 Smith-Waterman 算法 [ 27 ]。

RT-qPCR 验证

通过定量 PCR 验证了干旱响应 circRNA 的表达谱。使用Total RNA Kit (Tiangen, Beijing, China)根据制造商的说明从叶子中提取用作模板的总RNA，使用PrimeScriptTM RT试剂盒从20μl体积的1000ng总RNA合成第一条cDNA链根据制造商的协议，使用 gDNA 橡皮擦（Perfect Real Time）（Clontech，Shiga，Japan）。从差异表达的 circRNA 中随机选择 8 个 circRNA，并使用 RT-qPCR 进行分析。使用 Primer5 软件设计引物 [ 28]，并在 7500 实时 PCR 系统（Applied Biosystems，CA，USA）上进行 RT-qPCR。总反应体积为 20 μl，包含 10 μl 2X SYBR Premix Ex Taq?（TaKaRa Bio Inc.，日本）、1 μl 互补 DNA (cDNA) 反应混合物、0.5 μl 每种引物、0.5 μl ROX Reference DyeII 和 7.5 μl ddH 2 O. 3-磷酸甘油醛脱氢酶 (GAPDH) 基因用作标准化的看家基因 [ 29 ]。我们的 PCR 实验中使用的引物序列在S1 文件中进行了描述。PCR方法如下：95°C预变性30 s，95°C变性3 s，60°C退火30 s，55-95°C进行熔解曲线分析。所有反应均使用生物学一式三份进行。2 -ΔΔCT方法用于计算对照和处理植物之间基因表达的相对变化[ 30 ]。

结果

梨中circRNA的鉴定

为了研究参与干旱胁迫的 circRNA，我们进行了一项 RNA-seq 实验。从暴露于干旱胁迫的桦叶梨植物和对照植物的叶子构建了六个 cDNA 文库。在去除接头和引物序列以及短的低质量序列后，我们总共获得了 854,138,996 个干净的读数。Q20 和 Q30 得分均大于 90%（表格1)，说明测序质量高。经序列分析，得到889条平均长度为5,375 bp的非冗余circRNA（S2文件），其中614条（69.07%）的长度在150~5,000 bp之间，124条（13.95%）的长度为5,000 到 10,000 bp，所有其他长度都大于 10,000 bp (图。1）。Clean read 已存放在国家生物技术信息中心 (NCBI) 数据库中，可通过登录 SRP 150567 访问。

图1 梨形circRNA的长度分布。

表格1 通过 RNA-seq 获得的六个文库的测序数据。

在鉴定的circRNA中，186个（28.65%）是由外显子（exonic circRNAs）产生的（图2)，173 个 (24.26%) 是基因间的，70 个 (7.05%) 是由内含子产生的，其余的是外显子_内含子。

图2 四类型 circRNAS 的百分比。

circRNAs 响应干旱胁迫处理的差异表达

基因的表达模式可用作其推定生物学功能的指标。为了确定哪些 circRNA 在对照组和治疗组之间差异表达（每个条件三个生物学重复），根据统计阈值（FDR ≤ 0.05 和 〡log2 (ratio)〡≥ 1）过滤 circRNA。发现与对照组相比，治疗组中有 33 种 circRNA 的表达水平显着不同。10 个基因在干旱胁迫下下调，而 23 个基因上调（S3 文件）。差异表达的基因使用热图进行可视化（图 3）。

图 3 与对照相比，干旱胁迫下差异表达 circRNA 的热图。

参与干旱胁迫的circRNA的功能分析

为了进一步了解circRNA的潜在功能，我们预测并分析了它们的宿主基因。将 33 个差异表达的 circRNA 的宿主基因分为 171 个功能组（S4 文件），分为三个主要的 GO 分类类别（“生物过程”、“细胞成分”和“分子功能”），其中包含 100、18、和 53 个官能团，分别。氧化还原过程 (GO: 0055114) 由 11 个基因组成，在生物过程类别中占优势。具有 6 个基因的质膜 (GO: 0005886) 在细胞成分中占主导地位。Top GO富集分析结果显示在图 4.

图 4 差异表达的circRNA宿主基因的GO功能分析。

我们进一步分析了差异表达的 circRNA 的宿主基因以进行 KEGG 通路富集，并根据它们在 10 条通路中的参与情况对其进行注释（S5 文件）。在这些通路中，三个宿主基因被分配到“代谢过程”，而其他通路则与一两个基因相关。图5）。最常见的途径包括“次级代谢物的生物合成”、“核糖体”、“内质 中的蛋白质加工”、“泛素介导的蛋白水解”、“RNA 转运”。

图5 差异表达的circRNA宿主基因的KEGG分析。

circRNAs的靶miRNAs

最近的研究表明，circRNAs可以结合miRNAs，阻止它们靶向mRNAs，从而调节基因表达[ 13 ]。为了识别靶向 miRNA 的梨 circRNA，我们使用 psRobot 软件来预测潜在的 miRNA 结合位点。总共有 309 个 circRNA 被预测作为相应 180 个 miRNA 的海绵（S6 File）。在这 309 个 circRNA 中，146 个具有两个以上的 miRNA 结合位点。这意味着这些 circRNA 可以作为 miRNA 海绵。circRNA322有15个miRNA结合位点：miR5658、miR169a-5p、miRNA169b-5p、miR169c、miR169d、miR169e、miR169f-5p、miR169g-5p、miR169h、miR1616、miR1616、miR1619、miR1619、miR169i 这表明单个 circRNA 可以靶向多种 miRNA，并且单个 miRNA 可以被不同的 circRNA 靶向。例如，miR414 可以被 20 个 circRNA 靶向，而 miR-156 可以被 52 个 circRNA 靶向。33个差异表达的circRNA的潜在靶miRNA部分显示在图 6.

图 6 差异表达的 circRNA 及其目标 miRNA 的 络。

RT-qPCR 验证

为了验证转录组基因表达谱的可靠性，随机选择了 8 个差异表达的 circRNA 通过 RT-qPCR 进行表达分析。图 7）。RT-qPCR 结果中显示的表达模式与 RNA-seq 结果一致。例如，干旱胁迫后circRNA527的相对表达增加，而circRNA822的表达降低，这与RNA-seq结果一致。表明本次实验数据分析的结果是可靠的。

图 7 通过 RT-qPCR 测量的八种选定 circRNA 的相对表达。

讨论

过去一年，circRNA被认为是RNA剪接错误[ 31 ]。最近的研究已经在哺乳动物中鉴定出大量的 circRNA，并证明天然 circRNA 在动物的生物学和发育过程中发挥着重要作用 [ 10 , 32 ]。然而，与动物 circRNA 相比，植物 circRNA 很少受到关注 [ 17 , 33 ]。在这里，我们 告了对桦木梨中circRNAs的首次鉴定和表征，我们总共获得了889个circRNAs，平均长度为5,375 bp。在其他植物中，先前在拟南芥中鉴定了 6,012、12,037、854、113 和 496 个 circRNA [ 18]。]、水稻[ 17 ]、番茄 [ 14 ]、狗尾草斜体和玉米[ 34 ]。我们在确定梨的circRNAs数明显多于番茄，已经确定更大的小 。斜体和Z。mays，但比拟南芥和O少。苜蓿。在拟南芥和O . 苜蓿, cDNA 文库由多种组织类型生成；相反，本研究中的 cDNA 文库仅从叶子中生成。因此，我们可能只鉴定了桦叶梨中总 circRNA 的一个子集。

由于 circRNA 没有 5? 和 3? 末端，相关的内在稳定性使它们成为面对环境挑战时维持体内平衡的有力候选者 [ 35 ]。circRNA 也经常遵循组织和阶段特异性表达模式 [ 11 , 36 ]。在人类中，与主要的典型线性同种型相比，数据的所有循环转录本都以低水平表达 [ 37 ]。circRNA 子集的成员在年老小鼠与年轻小鼠的脑组织中显着上调，而一些则下调 [ 38 ]。circRNA 也已被证明在动物大脑 [ 39 ]、非酒精性脂肪性肝炎样本中高度丰富并动态表达 [ 39 ]。40 ]、维管细胞 [ 41 ]、各种胁迫下的拟南芥[ 42 ] 和各种水稻组织 [ 12 ]。在这项研究中，我们发现 33 个 circRNA 在脱水胁迫下差异表达，因此可能在梨的干旱胁迫耐受性中发挥重要作用。circRNA 的明显调控似乎是一种普遍现象。

植物对干旱胁迫的反应是一个复杂的过程，涉及多个基因和代谢 络，激素合成是重要因素之一[ 34 ]。2 宿主基因与应激反应相关。氧化还原加工途径是干旱胁迫的另一重要途径，本研究共有11个宿主基因与该途径相关。我们还鉴定了 1 个与泛素相关的宿主基因，这些基因可能参与信 转导和蛋白质在应激反应中的降解 [ 43 ]。“代谢过程”是与“生物过程”相关的 100 个亚组中的主要亚类。这一结果与之前的研究一致 [ 44 , 45]。这为进一步了解circRNA在干旱响应中的作用提供了进一步的见解。

CircRNA 可以作为竞争性内源 RNA 来结合 miRNA。在哺乳动物中，ciRS-7 有超过 70 个潜在的 miR-7 miRNA 结合位点。在番茄中，61 个 circRNA 充当了 47 个 miRNA 的海绵 [ 46 ]。在我们的研究中，预测 309 个 circRNA 充当 180 个 miRNA 的海绵。因此我们可以推断，包含常见 miRNA 结合位点的各种 circRNA 可能充当 miRNA 海绵来调节梨对干旱胁迫的反应。

总之，我们在桦叶梨中鉴定了 889 个 circRNA，其中 33 个 circRNA 显示出脱水响应。功能分析表明，差异表达的circRNAs参与了许多脱水反应过程，如代谢途径、内质 中的蛋白质加工和氨基酸的生物合成。circRNA-miRNA共表达 络表明circRNA参与了干旱响应过程。我们的研究结果为发现与干旱胁迫相关的基因提供了丰富的遗传资源，并且可以很容易地应用于其他果树物种。