1. 引言
Cas9是II型CRISPR/Cas (Clustered regularly interspaced short palindromic repeats and CRISPR associate)系统的核酸内切酶，可在单向导RNA (single-guide RNA, sgRNA)分子的引导下，特异性地识别靶标DNA，并对其进行双链切割 [1] - [4] 。由于其操作简便、切割效率高等优点，CRISPR-Cas9系统被广泛应用于基因组编辑、基因表达调控等领域 [5] - [7] 。近年来，国内外也已开展了一系列基于该系统的基因治疗基础研究，利用最为成熟的SpCas9系统进行白内障、地中海贫血、血友病等单基因疾病的基因纠正 [8] - [10] 。但是由于SpCas9编码序列长、蛋白过大，无论通过病毒载体进行基因转移，还是直接进行蛋白穿膜运载，都存在着较大的靶细胞运输难度，这限制了CRISPR/Cas9技术在基因治疗领域的应用 [11] [12] 。为解决这些问题，亟需在保证SpCas9活性与功能的同时，对其进行截短改造，拓展其应用前景。
国内外已通过结构与功能研究部分解析了SpCas9蛋白的结构及其切割机制。其中，SpCas9蛋白由1368氨基酸残基构成的单链折叠而成，主要组成两个叶片状(lobe)结构：REC-lobe主要负责靶DNA的识别，NUC-lobe负责靶DNA的剪切，两个lobe中间的凹槽为SpCas9、sgRNA和靶标DNA三者的相互作用提供空间。图1展示了已报道的SpCas9-sgRNA-DNA复合物三维结构模型(PDB ID: 4UN3) [15] ，可以
Figure 1. The atomic structure of the SpCas9-sgRNA-DNA ternary complex [15]
图1. SpCas9-sgRNA-DNA三元复合物原子结构 [15]
看出SpCas9的两个lobe还可以分为多个次级结构域，它们之间相互协调，共同实现对DNA的识别和剪切。具体来看：REC-lobe首先与sgRNA结合，通过sgRNA的碱基互补配对识别并结合靶标DNA；接着，SpCas9 NUC-lobe的两个核酸酶次级结构域(HNH及RuvC结构域)分别对互补及非互补DNA链进行剪切 [13] - [16] 。在这个过程中，SpCas9的构象也会发生变化。次级结构域的存在使SpCas9易于在结构域水平进行内部的局部优化，其改造难度降低。为了保证酶的完整活性，必须对结构域间的相互作用关系进行详细分析。但是，目前关于这些次级结构间作用关系的信息非常缺乏。
针对上述问题，本文选用分子动力学(Molecular dynamics, MD)模拟和动态网络分析方法研究SpCas9-sgRNA-DNA复合物结构，以三者间相互作用关系为基础分析SpCas9蛋白次级结构域间的相互作用，最后找出SpCas9蛋白次级结构域间相互作用的关键氨基酸 [17] [18] ，以加深人们对SpCas9作用机制的了解，也为后续的SpCas9蛋白的优化改造提供理论指导。
2. 方法和原理
2.1. SpCas9-sgRNA-DNA三元复合物长程分子动力学模拟分析
2.1.1. 模型构建
我们以目前完整程度最高的SpCas9复合物晶体结构(PDB ID: 4UN3)为基础，构建完整的复合物原子模型，用于后续MD模拟和动态网络分析 [15] 。由于4UN3中SpCas9蛋白和DNA链都有片段缺失，所以我们首先使用软件补齐缺失片段。对于SpCas9缺失的5个8~17残基长度的肽段，使用Modeller通过同源模建的方法进行补齐 [19] 。随后用AmberTools的NAB模块进行非互补DNA链的延长 [20] [21] 。由于4UN3结构中只包含11个碱基(nucleotide, nt)的非互补DNA链，且不含切割部位，因此我们新建了6 nt的B-DNA片段，与原非互补DNA链连接。新建DNA链穿过SpCas9切割活性位点。为保证新建DNA链处于天然态构象，通过NAMD程序对延长片段进行能量最小化处理 [22] 。SpCas9蛋白补齐片段和DNA链延长片段的序列信息均来自相关文献 [15] 。
将新建的全长SpCas9、非互补DNA链与原有sgRNA、互补DNA链组装，使用NAMD的IMD功能对新建片段进行局部微调，以消除新引入组分的构象冲突 [23] 。接着，在活性中心添加合适的催化金属离子(Mg2+)，获得完整的SpCas9-sgRNA-DNA复合物的结构原子模型。
2.1.2. 显式溶剂MD模拟
为研究复合物各个区域间的相互作用，我们用NAMD对新构建的复合物模型进行了长时间尺度的显式溶剂MD模拟。用分子可视化程序VMD建立模拟用的复合物—溶剂系统，采用CHARMM27力场描述复合物的蛋白质与核酸分子，溶剂水模型采用TIP3P模型 [24] [25] 。把复合物放置于一定大小的方形盒子中，保证其周围有足够数目的水分子，并加入一定数目的抗衡离子(Na+或Cl−)，以准确表示溶剂的离子浓度(0.15 mol/L l)及保持系统电中性。所构建的模拟体系盒子尺寸为160.8 Å × 151.2 Å × 186.5 Å，包含了约135,000个水分子。MD模拟基于恒温、恒压和恒定粒子数系综(N-P-T系综)和周期性边界条件进行：温度和压力分别为298.15 K和1 atm，模拟所用的时间步长为2 f。在模拟中，应用Berendsen方法调整系统的温度与压力，非键静电力采用PME方法计算，范德华力计算采用截断值方法(cutoff = 10 Å)，与氢原子相连的化合键的键长用SHAKE算法来约束。
2.2. 动态网络分析
2.2.1. 动态网络简介
动态网络是一种分析复杂网络的方法。近年来，随着小世界理论和无标度理论的发展，该方法已经逐渐被应用到生物学领域，并取得实质性进展，在结构生物学方面收获丰硕。该方法的中心思想是把蛋白质等生物大分子结构中的每个残基作为一个节点，通过特殊的数学算法，分析节点间的相互联系，相互作用的节点间以边连接，得出所有节点和边组成的动态网络图。分析此动态网络图隐含的信息，并推断这些信息在分子结构中对应的生物学信息 [17] [18] 。
在构建动态网络的过程中，首先根据模拟轨迹计算所有节点两两之间的相关系数 [26] 。相关系数的取值范围为−1~1，值越大说明相关性越强。随后找出作用比较密切的节点对，以边连接，并根据两个节点之间的相关系数计算边的强度。其中，密切作用节点对的判定标准为：如果两个节点间的距离在75%的模拟轨迹中都小于4.5 Å，则认为两个节点间作用密切。
动态网络图生成以后，进行社区结构划分。所谓社区，是指动态网络中的一个亚网络，在此亚网络中，内部节点间的连接稠密程度高于与外部节点的连接。事实上，社区之间也存在边进行连接，而这些边的存在对社区间的通讯起着至关重要的作用。与分子结构进行对照，就可以找出紧密联系的残基，划为独立的结构域，并找出结构域间相互作用的关键位点。
2.2.2. 复合物关键位点的动态网络分析
使用VMD的动态网络分析模块对SpCas9复合物模拟轨迹进行分析。首先以每个氨基酸和核苷酸为节点，构建所有节点两两之间相关系数的分布图；基于此相关系数分布图，计算动态网络图，进行社区的划分；最后，找出所有社区间的关键节点，把SpCas9蛋白分别与sgRNA、DNA作用的关键节点进行RMSD分析，验证关键位点可靠性。
对复合物模拟轨迹进行RMSD分析的理论依据为：在分子动力学模拟轨迹中，强烈相互作用的残基在系统稳定后会始终维持比较固定的相对构象，那么这些残基组合的RMSD值通常保持在较小的水平。具体思路如下：根据动态网络分析的原理，相互作用关键节点的距离在75%的模拟轨迹中都小于4.5 Å，意味着这两个节点代表的残基始终保持稳定的相互作用关系，那么此残基对的RMSD值将维持在较小水平。例如，分别以a、b表示这两个残基，其中a为氨基酸残基，c表示a前后四个残基中与a结构类似，且在网络分析中与b不存在连接边的残基。如果动态网络分析准确，那么a、b间作用会比c、b间作用稳定，则ab组合的RMSD值小于cb组合的RMSD值。
3. 结果与讨论
3.1. 三元复合物原子模型构建与MD模拟分析
三元复合物原子模型构建过程及新建复合物结构如图2所示。其中，补齐的Cas9缺失区域未对Cas9蛋白整体结构造成影响，而延长的非互补DNA链则包围在Cas9切割活性中心区域。经进一步与sgRNA、互补DNA链组装后，获得了完整的复合物模型。
我们对上述复合物进行约100 ns的MD模拟。然后，从每一帧MD模拟轨迹中抽提出复合物的构象，并以第一帧为参考构象计算后续每一构象的RMSD值，结果如图3所示。由图可见，整个体系在约30 ns后趋于稳定。因此，后续的动态网络分析使用的是30 ns以后的MD模拟轨迹。
3.2. MD模拟轨迹的动态网络分析
MD模拟轨迹中三元复合物中SpCas9残基、sgRNA及DNA的相关性如图4所示。其横坐标和纵坐标均表示复合物的各个残基单元，每个点的颜色由对应残基间的相关系数决定，颜色越浅，值越大。其
Red dash in SpCas9 lines are missing fragments in the crystal structure. Residues shown as dots are active sites for shearing non-complementary DNA strand. 红色虚线表示SpCas9晶体结构中缺失片段。点状残基表示切割非互补DNA链的活性位点。
Figure 2. Schematic diagram of building atomic model of SpCas9-sgRNA-DNA complex
图2. SpCas9-sgRNA-DNA复合物结构原子模型的构建过程
Figure 3. MD RMSD of the ternary complex with respect to the initial structure
图3. MD模拟轨迹的复合物RMSD值
Figure 4. MD correlation coefficients of the ternary complex
图4. 复合物模拟轨迹的相关系数热图
中，左下–右上方向的浅色对角线表示每个残基与其自身及邻近残基存在高度相关性。以这条对角线为轴，相关系数对称排列，与残基间相关系数的计算顺序相符。图中的相关系数分布明显分为两个区域，分别表示核酸结构与蛋白结构。其中核酸区域又分为两部分，分别为sgRNA和DNA，蛋白区域也分为两部分，分别为SpCas9的REC-lobe和NUC-lobe，这些分布都与复合物的结构高度吻合，也说明相关性分析的可靠性。
图5展示了基于相关系数分布进行动态网络分析的结果。图5(a)是所有节点间的动态网络图，相关节点间以边连接。此网络图与复合物结构吻合。进一步地，将节点间的边根据作用强弱进行加权(图5(b))，可以看出，同一个结构域内的节点间作用更为强烈。我们对加权的网络进行了社区划分，结果如图5(c)所示，各个社区被标记不同颜色。与复合物结构对照可以发现，社区的划分与三元复合物中sgRNA、DNA及SpCas9的次级结构域高度吻合，因此通过分析社区间的相互作用，可以揭示三元复合物各结构域的相互作用。为找出SpCas9结构域间相互作用的关键氨基酸，我们进行了社区间相互联系的关键节点的计算(图5(d))，并将其叠合到三元复合物原子结构模型中。可以发现，多数节点作用对都位于结构域之间，符合次级结构域相互作用的空间位置。
我们发现，在所有相互联系的关键节点中，共包含7组氨基酸-核苷酸作用节点。SpCas9晶体结构也显示这7个氨基酸与sgRNA/DNA发生相互作用 [14] 。为验证关键节点分析可靠性，我们对上述7个氨基酸–核苷酸对进行了RMSD分析。表1列出了动态网络分析获得的氨基酸-核苷酸节点对，及RMSD分析过程中用到的对照氨基酸。对照氨基酸均位于关键氨基酸上下两个位点以内，且结构性质与关键氨基酸相似。MD稳定轨迹中各关键氨基酸及对照氨基酸全原子与相应核苷酸的RMSD值及其统计结果如图6所示。在RMSD计算过程中以MD模拟稳定后的第一帧轨迹作为参考构象。
其中，五组节点对的RMSD值显著低于对照组，表明节点对间存在强烈的相互作用。仅有49ADE-76LYS、45URA-401LYS两个节点对的RMSD值稍大，其原因可能为测试节点氨基酸(赖氨酸)与对照组(精氨酸)的侧链转动能力不同。分析MD模拟轨迹发现，尽管赖氨酸侧链末端的氨基与相互作用的脱氧核糖核酸在绝大多数时间内保持稳定构象，但由于长侧链活动性很强，末端氨基少量偏移就可使整个侧链发生很大构象变化。而精氨酸侧链只包含三元甲基链，且在侧链末端是一个体积巨大的胍基，胍基的存在使甲基链几乎没有空间进行转动，导致精氨酸侧链的构象在整个MD模拟过程中很少变动。所以，赖氨酸与精氨酸相似的RMSD中值已以说明赖氨酸与脱氧核糖核酸作用的稳定性，这也表明7对氨基酸–核苷酸在模拟中始终保持强烈稳定的作用关系。此外，网络分析得到的832 Arg和1311 His也被实验证实与非互补DNA链的识别效率有关 [27] 。这些结果也说明，通过动态网络分析识别关键相互作用节点对有相当的可靠性。
基于动态网络分析结果，我们将得到的关键互作位点进行了分类，结果如表2、图7所示。其中，表2列出了序列间隔大于10的氨基酸-氨基酸关键节点对。这些节点对可以分为三个部分：第一部分中各节点对包含的两个节点均位于NUC-lobe，为切割结构域内部的关键氨基酸，主要与DNA切割有关；第二部分中各节点对包含的的两个节点均位于REC-lobe，为识别结构域内部的关键氨基酸，主要与DNA识别有关；第三部分中各节点对包含的两个节点分别位于不同lobe，在识别到切割的转换过程中起重要作用。计算分析结果表明，这些位点为结构域间联系的关键位点，十分有必要在蛋白截短优化中保留。这些位点的详细功能还有待进一步的实验验证。
4. 结论
为揭示SpCas9各结构域间的相互作用，本文在已有晶体结构基础上，构建出完整的、切割活性状态下的SpCas9-sgRNA-DNA三元复合物结构；通过MD模拟和动态网络分析，进行复合物各残基间相互联
Figure 5. Dynamical network analysis. (a) Unweighted network of SpCas9 complex; (b) Network of SpCas9 complex shown with edge widths corresponding to their weights; (c) Communities of SpCas9 complex painted with different colors; (d) Key node pairs connecting interacting communities. Key nodes are shown as yellow dots. Interacting amino acid-amino acid nodes and amino acid-nucleotide nodes are connected with yellow and red lines respectively
图5. 动态网络分析。(a) SpCas9复合物无权重网络；(b)边宽与权重相对应的SpCas9复合物网络；(c)以不同颜色表示的SpCas9复合物社区划分；(d)连接互作社区的关键节点对。关键节点以黄点表示。互作氨基酸–氨基酸节点与互作氨基酸–核苷酸节点分别以黄线和红线连接
Table 1. Amino acid-nucleotide key node pairs
表1. 氨基酸–核苷酸关键节点对
RMSD of key amino acid-nucleotide is colored in cyan while RMSD of key amino acid-nucleotide in control group is colored in pink. 实验组与对照组关键氨基酸-核苷酸的RMSD曲线分别以蓝绿色和粉红色表示。
Figure 6. RMSD analysis of amino acid-nucleotide key node pairs
图6. 氨基酸–核苷酸关键节点对RMSD分析
Table 2. Key pairs of SpCas9 amino acids
表2. SpCas9氨基酸关键节点对
Each amino acid in a pair is labeled in the same color. 每组关键节点对的两个氨基酸以相同颜色表示。
Figure 7. Schematic of key pairs of SpCas9 amino acids
图7. SpCas9氨基酸关键节点对示意图
系的动态网络分析及社区划分，发现网络社区与复合物三个组分及SpCas9次级结构域间高度吻合；最后，通过社间相互联系的关键节点的分析，获得了相互作用的关键氨基酸。通过RMSD分析及与已报道实验结果对照，所预测的关键氨基酸信息可靠。因此，本研究不仅可以加深人们对SpCas9作用机制的了解，而且为后续的SpCas9的优化设计提供重要的理论指导。
致谢
感谢复旦大学生命科学学院袁慧同学对文本修改提出宝贵意见。本工作的计算研究得到国家自然科学基金(91430112)、NSFC-广东联合基金(第二期)超级计算科学应用研究专项资助和国家超级计算广州中心支持。
*通讯作者。