key: cord-0075979-toj8det3 authors: Zhuang, Haiming; Fan, Xiaohua; Ji, Dong; Wang, Yuanhao; Fan, Jigang; Li, Mingyu; Ni, Duan; Lu, Shaoyong; Li, Xiaolong; Chai, Zongtao title: Elucidation of the conformational dynamics and assembly of Argonaute–RNA complexes by distinct yet coordinated actions of the supplementary microRNA date: 2022-03-07 journal: Comput Struct Biotechnol J DOI: 10.1016/j.csbj.2022.03.001 sha: 12655a48d0ae836b299c69b72b754c60623c181a doc_id: 75979 cord_uid: toj8det3 Argonaute (AGO) proteins, the core of RNA-induced silencing complex, are guided by microRNAs (miRNAs) to recognize target RNA for repression. The miRNA–target RNA recognition forms initially through pairing at the seed region while the additional supplementary pairing can enhance target recognition and compensate for seed mismatch. The extension of miRNA lengths can strengthen the target affinity when pairing both in the seed and supplementary regions. However, the mechanism underlying the effect of the supplementary pairing on the conformational dynamics and the assembly of AGO–RNA complex remains poorly understood. To address this, we performed large-scale molecular dynamics simulations of AGO–RNA complexes with different pairing patterns and miRNA lengths. The results reveal that the additional supplementary pairing can not only strengthen the interaction between miRNA and target RNA, but also induce the increased plasticity of the PAZ domain and enhance the domain connectivity among the PAZ, PIWI, N domains of the AGO protein. The strong community network between these domains tightens the mouth of the supplementary chamber of AGO protein, which prevents the escape of target RNA from the complex and shields it from solvent water attack. Importantly, the inner stronger matching pairs between the miRNA and target RNA can compensate for weaker mismatches at the edge of supplementary region. These findings provide guidance for the design of miRNA mimics and anti-miRNAs for both clinical and experimental use and open the way for further engineering of AGO proteins as a new tool in the field of gene regulation. Argonaute (AGO) proteins play a critical role in all small-RNAguided RNA silencing, pivotal in the regulation of gene expression and essential for cell proliferation, differentiation, and homeostasis [1] [2] [3] . MicroRNAs (miRNA), which are short non-coding RNAs of $ 22 nucleotides, function by guiding AGO proteins to interact with complementary target RNAs, yielding the RNA-induced silencing complex (RISC) [4, 5] . The formation of RISC can promote the translation repression of messenger RNA [6, 7] or the degradation of exogenous pathogen RNA [8] . In addition to the key role in cellular process, microRNA has also been widely used as an experimental tool [9, 10] and recent advances have gradually made miRNA therapeutics for clinical utility [11] [12] [13] [14] . There are four kinds of AGO proteins in human cells: AGO1, AGO2, AGO3, and AGO4. Among them, AGO2 is known as the only slicer [15] . AGO2 is a bi-lobed protein with a central cleft for binding of guide and target RNAs and is composed of four globular domains (N, PAZ, MID and PIWI) and two connective linker domains (L1 and L2) [16] (Fig. 1A and B) . The 5 0 and 3 0 ends of miRNA are anchored at the MID and PAZ domain, respectively [17] . In the crystal structure of miRNA and target RNA bound to the AGO2 protein, the guide miRNA can be divided into five func-tional domains, including anchor (g1, from the miRNA 5 0 end), seed (g2-g7 or g2-g8), central (g8-g13 or g9-g13), 3 0 supplementary (g13-g17), and tail (g18-the end) [18] (Fig. 1C-E) . In miRNA targeting, initial pairing is formed in the seed region [19] . Seed region pairing alone is sufficient for target recognition [20] and $ 80% of miRNA interactions rely on seed complementarity [21] . Seed region is also the most evolutionarily conserved domain in animals [22] . Owing to the evolutionary conservation of seed region, miRNA family members share identical seed sequences. However, miRNAs in the same family have different targets derived from the determinant role of 3 0 end pairing in target specificity [23, 24] . The 3 0 end of miRNA can be divided into 3 0 supplementary and tail [18] (Fig. 1D) , and the supplementary pairing beyond seed region could efficiently regulate target mRNA levels [25] . The additional supplementary interaction can enhance target recognition or compensate for seed mismatch [26, 27] . Previous studies showed that the additional supplementary pairing could enhance target RNA binding affinities with varied abilities dependent on different miRNAs: $2-fold for let-7 [18] , $7-fold for miR-21 [28] and >20fold for miR-122 [29] . Biochemical experiments and mathematical model further elucidated that AGO proteins were highly occupied by seed-plus-supplementary-pairing targets even with the existence of overnumbered seed-pairing-only targets, revealing the differential targeting driven by supplementary pairing [27] . For the biochemical properties of supplementary pairing, catalytic kinetic analysis showed that two simultaneous mismatches (g13 and g17) introduced into pairing slightly weaken the target cleavage efficiency [20] . Competition assays also confirmed that g15 and g16 mismatches did not impair binding affinity of target RNA in mouse AGO2-RISC [18] . Consistent with the experimental results, target predicting algorithms implied that 3 0 pairing effect was relatively insensitive to edge mismatches while sensitive to continuous Watson-Crick pairing [23] . In addition to the supplementary region, the miRNA's 3 0 tail region can pair with the target RNAs. This pairing paradigm in the tail region based on the seed and supplementary pairings could lead to the 3 0 end destabilization of miR-NAs because of the release of 3 0 end from the PAZ domain [30, 31] . When in vivo instability of miRNAs, the guide miRNA may decay instead of target repression; this biological phenomenon is referred to as target-directed miRNA degradation (TDMD) [32] [33] [34] . MiRNA can exhibit heterogeneous 3 0 ends due to the posttranscriptional 3 0 addition of nucleotides which may stabilize miRNA or strengthen guide-target interaction [35, 36] . Recent biochemical data showed that for the seed-plus-supplementary pairing, the increased lengths of miRNA from 21mer to 23mer could enhance target affinity [29] . In contrast to seed-plussupplementary pairing, only-seed pairing had no appreciable difference on target RNA affinities with extended lengths of miRNA. Based on the crystal structure of seed-plus-supplementary pairing AGO2 complex, the structural model with the supplementary pairing shows that the PAZ domain and the L2 loop shift away from the central cleft, thereby providing sufficient space to form a supplementary chamber for guide-target pairing [29] . The movement of the PAZ domain may create tension in the miRNA 3 0 tail, and the increased target RNA affinity may be ascribed to the relieved tension at the miRNA 3 0 tail, which facilitates the opening of supplementary chamber. However, the underlying mechanism how the structure of the supplementary chamber and the AGO-target interactions influenced by the miRNA extensions remains unclear. Furthermore, the structural model of the supplementary pairing is unable to elucidate the mechanism how the mismatches at the edge of supplementary pairing had a minor effect on target binding and the 3 0 pairing was relatively insensitive to the predicted thermodynamic stability. Although the crystal structure of supplementary pairing in the AGO complex has been solved, it is challenging to provide enough information for elucidating the mechanism of supplementary pairing with only one averaged snapshot because of the limited conformational changes. Molecular dynamics (MD) simulations that explore protein and nucleic acid conformational dynamics at the atomic level are critical for understanding of target repression and TDMD process [37] [38] [39] [40] [41] [42] [43] . In the recent work, we have used MD simulations to directly uncover biomolecular mechanisms and protein-ligand/protein or protein-DNA/RNA recognitions [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] . Previous investigations on AGO utilizing large-scale MD simulations have successfully illuminated detailed mechanisms behind AGO-RNA interaction [58] [59] [60] and have been testified by experiments [61] . Thus, this method can further help to improve miRNA target site prediction algorithm which can improve the identification of genes directly regulated by miRNAs [62, 63] and to guide rational engineering of AGO complex for therapeutical implications. Here, we performed a multiple, microsecond-scale length explicit MD simulations to investigate the effect of different miRNA lengths and miRNA-target pairing patterns on the target RNA binding affinities. Nine systems were built with different miRNA lengths and miRNA-target pairing patterns. The simulations revealed that with the extension of the miRNA lengths, the supplementary chamber opened a smaller range and fastened the mouth of the central cleft to prevent the escape of target RNA from the complex and protect it from the attacking from solvent. Our study further showed that the inner stronger matching pairs between the miRNA and target RNA can compensate for weaker mismatches at the edge of supplementary region and result in almost the same pairing effect between the guide and target RNAs. This study aims to shed light on the mechanistic understanding of the supplementary pairing paradigms, inspiring further research on target prediction algorithm and miRNA design for experimental and therapeutic use. A total of nine systems were performed: three for seed-pairing (g2-g8)-only (T1 systems), three for seed-plus-supplementary (g13-g16)-pairing (T2 systems), and three for seed-plussupplementary (g15-g19)-pairing (T3 systems). For systems with the same pairing pattern, miRNA length was 21mer, 22mer, and 23mer, respectively. The RNA sequences of the nine systems were shown in Fig. 1C -E. The crystal structure of seed (g2-g8)-plus-supplementary (g13-g16)-pairing AGO2 complex (PDB ID 6N4O) [29] was selected as the initial structure for further modelling. The missing residues and atoms of AGO protein were added using the MODELLER program except the first 21 residues at the N-terminus. The missing g10 nucleotide of the guide RNA was filled up using the counterpart of the guide RNA based on the seed-pairing-only structure of AGO2 complex (PDB ID 4Z4D) [64] because of the similar central region of guide RNA in the two structures. The missing g19 was added through growing a nucleotide from g18 and ligating it with g20 using the Discovery Studio 2019. The missing t8 nucleotide in the target RNA was added by the counterpart of the target RNA from the TDMD structure of AGO2 complex (PDB ID 6NIT) [30] . After that, a nucleotide was prepended from t8 to t7 and t6 was added by lengthening from t5 using the Discovery Studio 2019. To remove the incompatibility of newly added nucleotides, 10,000 steps of minimization of nucleotides were performed using the steepest descent algorithm with the protein fixed. Then, a second-round minimization was performed without any restraints. Three seed-pairing-only systems, including 21mer-T1, 22mer-T1, and 23mer-T1, were constructed. Based on the initial structure, mutations of nucleotides to form seed-pairing-only nucleic sequences and extensions of guide RNA from 21mer to 23mer were done using the Discovery Studio 2019. 10,000 steps of minimization were then performed using the steepest descent algorithm to optimize complex structures. Three seed-plus-supplementary (g13-g17)-pairing systems, including the 21mer-T2, 22mer-T2, and 23mer-T2, were constructed. First, to model the 21mer-T2 system, nucleotide mutation, nucleotide extension, and energy minimization were performed based on the initial structure. For the additional pairing of g17-t6 compared, a 20 ns MD simulation was performed to drive the guide and target pairing at the anticipative position. The prepared complex was solvated to a water box and Na + and ClÀ were added to simulate normal saline environment. The system was then minimized and heated. In the production simulation, a pulling force was added of 3 kcal/mol/ Å 2 to promote the pairing of g17 to t6 by restraining the distance between the N1 atom of g17 and the N3 atom of t6 and the distance between the N6 atom of g17 and the O4 atom of t6 within 3.5 Å. Next, to choose the most optimal conformation, average-linkage algorithm was carried out to extract the most representative structure based on the snapshots from the last 5 ns of simulations. We selected the coordinates of Ca atoms in the AGO2 protein and P atoms in the nucleotides as input for average-linkage algorithm. With the obtained structure of the 21mer-T1, the guide extension from 21mer to 23mer and energy minimizations were accomplished using the Discovery Studio 2019. Three seed-plus-supplementary (g15-g19)-pairing systems, including the 21mer-T3, 22mer-T3, and 23mer-T3, were constructed. First, to model the 21mer-T3 system, we mutated and remove redundant nucleotides using Discovery Studio 2019 from the structure of 21mer-T2. To enable the additional pairing of g18 and g19 in the 21mer-T3 compared to the 21mer-T2, a 20 ns simulation was performed as mentioned above. Pulling force was added using 3 kcal/mol/Å 2 to promote the pairing. After a 20 ns simulation, the 3 0 end of guide RNA was released from the PAZ domain. This phenomenon might due to the inappropriate force added. Since the pairing in the tail region may cause the release 3 0 end of guide RNA, the 21mer RNA was too short to bear the tension upon additional 3 0 pairing. Therefore, we chose 23mer-T2 to construct 23mer-T3 using the above method. We used averagelinkage algorithm to extract the representative structure as the initial structure of the 23mer-T3 system. With the obtained structure of the 21mer-T3, the guide truncation from 23mer to 21mer and energy minimization were accomplished using the Discovery Studio 2019. The above-mentioned systems were prepared with the LEaP program using ff14SB and ff99OL3 to describe the RNA-protein complexes [65] [66] [67] . The prepared complex was solvated to a truncated octahedron transferable intermolecular potential three point (TIP3P) water box [68] and Na + and ClÀ were added to neutralize and simulate normal saline environment. Two rounds of minimization were performed in which the first round contained 5000 steps maximum minimization cycles with the protein and nucleic fixed and the second round contained 10,000 steps maximum minimization cycles with no constraints. Subsequently, all systems were heated from 0 to 300 K within 300 ps, followed by 700 ps equilibration running in a canonical ensemble (NVT). After all these preparations, 10 replicas of independent 1 ls simulations were performed with random velocities under isothermal isobaric (NPT) conditions, producing a total of 90 ls simulation trajectories for analysis. For the long-range electrostatic interactions, the particle mesh Ewald method [69] was used, while a 10 Å non-bonded cut-off was introduced for the short-range electrostatics and van der Waals interactions. Covalent bonds involving hydrogens were restrained using the SHAKE method [70] . Principal component analysis (PCA) widely used in describing the kinetic process during simulation is a technique that can transform a series of potentially coordinated observations into orthogonal vectors. Among these vectors, the first two principal component (PC1 and PC2) provide the dominant motions during simulation [71] . In this analysis, PCs were generated based on coordinate covariance matrix of Ca atoms in the AGO protein using every 10 frames in all nine systems and these collected frames were all projected on the first and second principal components. Generalized correlation (GC ij ) analysis was applied to calculate the correlated motion. Compared to traditional Pearson correlation analysis, GC ij analysis is independent of mutual orientations of the atomic fluctuations and enable the separation between linear and nonlinear contributions [72] . To reflect how much information of one atom's position is provided by another, Mutual Information (MI) was introduced and calculated by:. The right side of equation can be related to the more widely known measure of entropy which calculated by:. To calculate based on the correlation between pairs of atoms:. was further related to a more intuitive Pearson-like correlation coefficient GC ij which can be calculated by:. where d represents the dimensionality of x i and x j , which is 3 in our study. GC ij calculation was done by g_correlation tool in Gromacs 3.3 [73] with the coordinates of Ca in each residue as input. The inter-domain correlation was calculated by accumulation of GC ij between each residue in respective domains with a minimum threshold value of 0.50. Then, the accumulated value was normalized according to the maximum and minimum value in the counterpart in all nine systems. Dynamic network analysis was performed to reflect the motion connection using the Network View plugin in VMD [74] . In our analysis, the Ca atoms in the AGO2 protein, the N1 atoms in the uracil and cytosine, and the N9 atoms in the adenine and guanine were chosen as nodes to represent their corresponding residues. Edges were drawn between nodes whose distances are within a cutoff of 4.5 Å for at least 75% of simulation time. The edge between nodes was calculated by:. where i and j represent the two nodes. In addition, community, which means a combination of residues whose connections are stronger was calculated using the Girvan-Newman algorithm [75] through the gcommunities program in VMD. Markov state model (MSM) is a mathematical framework to describe the dynamics of time-series data [76] [77] [78] . It has been used as an invaluable tool in understanding the dynamics during MD simulation [79, 80] . In this analysis, MSM was used to differentiate conformational states during simulation and help to extract representative structures from each state. All MSM calculations were performed using the PyEMMA software. MSM transition matrix was first calculated based on the probability of transition between different states. Implied timescale (ITS) test was performed to check the Markovian property and choose the proper lag time(s). ITS as a function of s can be calculated by:. where k i means the eigenvalue gotten from MSM transition matrix of the i th process. Based on this, implied timescale plot can be drawn in which each curve represents an average transition time in one process. When the curve became approximately constant, the corresponding lag time is appropriate for following analysis and the system is memoryless (Markovian property). In this study, the lag time was set as 10 ns. Every frame of each system was projected on the free energy plot according to the characteristic vectors (CV). Then, K-centers algorithm was applied to cluster the two-dimensional conformations into 300 microstates in each system. Based on these 300 microstates, PCCA + algorithm was performed to divide the microstates into certain states (can be thought as metastates). With the divided states, the Chapman-Kolmogorov test was computed to further testify the property of MSM. Every microstate was labeled as a certain metastate. The frames near the K-centers (microstates) belonging to the certain state in energy landscape were extracted to constitute a new trajectory to represent the state. The structure that has the smallest cas (RMSD) with other frames in the trajectory was chosen as the representative structure. A total of nine systems were performed: three for seed-pairing (g2-g8)-only (T1 systems), three for seed-plus-supplementary (g13-g16)-pairing (T2 systems), and three for seed-plussupplementary (g15-g19)-pairing (T3 systems). For systems with the same pairing pattern, miRNA length was 21mer, 22mer, and 23mer, respectively. The RNA sequences of the nine systems were shown in Fig. 1C -E. Each system had 1 ls  10 independent runs with random initial velocities, leading to a cumulative simulation timescale of 90 ls. Such an extensive timescale has been proven efficient for exploring the biomacromolecule conformational landscapes and recognition processes. To explore the effect of RNA lengths on overall conformational dynamics, the Ca atoms rootmean-square deviation (RMSD) of the AGO2 protein relative to the crystal structure was calculated. No significant differences were found except for the 21mer-T3 system (Supplementary Fig. 1A-C) . Compared to the 22mer-T3 and 23mer-T3 systems, a $ 0.5 Å RMSD increase in the 21mer-T1 system was observed, indicating a conformational variation compared with the initial structure. To explore the conformational variation, we extracted the representative structure of the 21mer-T3 system from the 10 replicas of MD simulations and found the release of the 3 0 end of the miRNA from the PAZ domain in each replica ( Supplementary Fig. 1D ). This observation indicated that the length of the 21mer miRNA was insufficient to tolerate the increased tension created by the 3 0 tail pairing. To further uncover the influence of miRNA lengths on local conformational dynamics, the root-mean square fluctuation (RMSF) of each residue (represented by the Ca atoms for AGO2 and P atoms for miRNA) was calculated ( Supplementary Figs. 2 and 3) . To better show the differences, we subtracted RMSF of the 22mer and 23mer systems from the counterpart of 21mer systems in the three systems and projected the difference values on the protein structure ( Fig. 2A-I) . Inspection of the RMSF plot showed that the PAZ domain displayed the highest fluctuation and showed an increased plasticity with the extension of miRNAs, indicating the less constrained conformation of the AGO2 protein in response to the extended miRNAs. However, it should be noted that with the extension of the miRNA lengths from 22mer to 23mer, the RMSF of the PAZ domain decreased significantly in the region proximal to the miRNA tail. This can be attributed to the enhanced stability created by the additional binding of 23th nucleotide. In general, the RMSF variations in the PAZ domain induced by different miRNA lengths are originated from the combined actions of the tension release and the additional binding. To characterize the global conformational transition and essential dynamics in the simulations, principal component analysis (PCA) was performed. We projected all MD snapshots from 10 independent runs of each system on the two-dimensional plot according to the first two principal components (PC1 and PC2) ( Fig. 3A -I) and the domain motions along the PC1 (Fig. 3J ) and PC2 (Fig. 3K) were shown on the AGO2 protein. The pairing pattern can be distinguished based on the PC2 with the seed-pairing-only systems showing higher PC2 value. Along the PC2, the PAZ domain underwent large-amplitude motions towards the N domain and moved away from the L2 domain, which rendered the supplementary chamber changing from the ''open" to ''closed" conformations (Fig. 3K ). This observation supported the notion that the PAZ domain shifted away from the N domain to provide sufficient space to form supplementary chamber. By superimposing the averaged structures of the 23mer-T1 system which represented the seed-pairing-only conformation and the 23mer-T2 system which represented the seed-plussupplementary-pairing conformation, the shift of the PAZ domain concerted with the motion of the PC2 (Fig. 3L) . However, PC2 values had no variation with the extended miRNA lengths, suggesting that miRNA lengths cannot affect the ''open" to ''closed" motions. To further elucidate this observation, we calculated the distance between the E261 Ca atom from the PAZ domain and the D358 Ca atom from the L2 domain, representing the shifted range of the PAZ domain from the L2 domain ( Fig. 3L and M) . The calculated distances in the seed-pairing-only systems were relatively larger than those in the seed-plus-supplementary pairing systems. However, no significant differences were found with the extension of the miRNA lengths in the additional supplementary pairing systems, indicating that the extended miRNA lengths did not affect the ''open" to ''closed" motions of the AGO2 protein in the seedplus-supplementary pairing systems. Compared to the PC2, the PC1 values changed with the extended miRNA lengths. In the seed-plus-supplementary-pairing systems, PC1 values increased with the extended miRNA lengths from 21mer to 23mer while this property was not observed in the seed-pairing-only systems (Fig. 3A-I) . Along the PC1, both the PAZ and PIWI domains moved towards each other and yielded the supplementary chamber to open a smaller range, which can be called as the ''departing" to ''approaching" motions (Fig. 3J) . To further probe the influence of miRNA lengths on the ''departing" to ''approaching" motions, we projected all MD trajectories onto the two-dimensional surface according to the distances from the I61 Ca atom to the G836 Ca atom (d 1 ) and the E333 Ca atom to the center of mass of Ca atoms in the A603, G604 and D605 (d 2 ) (Fig. 4A-I and Supplementary Fig. 4 ). D 1 represents the distance between the PAZ and PIWI domain which correlates to the ''departing" to ''approaching" motions while d 2 represents the distance between the PAZ and N domain. Since the PAZ, PIWI, and N domains form the mouth of the supplementary chamber, d 1 and d 2 can reflect the supplementary chamber-opening degree. To differentiate the conformational states in the energy landscape, Markov state model (MSM) was applied (see Material and Methods) and only the states in the energy basin of free-energy landscape were analyzed. The implied timescale test and the Chapman-Kolmogorov test were calculated to testify the Markovian property ( Supplementary Figs. 5 and 6 ). The proportion, the representative structure, and the representative trajectory of each state were calculated based on the MSM. In addition to the free-energy diagram, violin plots of d 1 and d 2 were also drawn to better compare between different systems ( Supplementary Figures 7 and 8) . Compared to the seed-plus-supplementary pairing systems, the range of both d 1 and d 2 parameters were generally larger in the seedpairing-only systems, indicating the instability of the contacts between the PAZ, PIWI and N domains when pairing only in the seed region. The higher d 2 parameter in the seed-pairing-only systems also reflected the larger distance between the PAZ and PIWI domains, which further suggested a contraction of the PAZ and PIWI domains with the additional supplementary pairing. According to the states extracted using MSM in the seed-pairing-only systems, with the increase of the d 1 parameter, three states were observed: CX, CY, and CZ ( Fig. 4A-C) . The proportion and location of the three states did not change with the extended miRNA lengths and the projected area on the free-energy landscape was widely distributed and very similar in three systems, which reflected the conformational similarity of the AGO2 protein with different miRNA lengths in the seed-pairing-only systems. In the seed-plus-supplementary (g13-g17) pairing systems, two different states were observed: one with a smaller d 1 parameter means a closer contact between the N and PIWI domains (CL), and the other was opposite (CR) (Fig. 4D-F) . In the seed-plussupplementary (g15-g19) pairing systems, three different states were observed (Fig. 4G-I) . One had larger d 1 and d 2 parameters (CU), meaning a distant relationship between the PAZ, PIWI and N domains which only existed in the 21mer-T1 system. The other two states shared a similar d 1 parameter and from the high to low d 2 parameter was CM and CD, respectively, which merged into one state in the 22mer-T3 system. Except for the 21mer-T3 system whose 3 0 end of miRNA was released from the PAZ domain, both d 1 and d 2 parameters broadly decreased with the extended miRNA lengths in the seed-plus-supplementary-pairing systems, indicat- ing the shrinkage of the mouth of supplementary chamber encompassed by the PAZ, PIWI and N domains. By superimposing the representative structures of their respective dominant state in the T2 systems (CR for 21mer-T2 and 22mer-T2, and CL for 23mer-T2), we found that with the extended miRNA lengths, the PAZ, PIWI and N domains were getting closer to each other and wrapping around the target more tightly (Supplementary Figure 9 ). To show the conformational changes clearly, we superimposed the representative structures of the CL, CR, CM, and CU states. Compared to the CR conformation, the proximity of the PIWI and N domains in the CL conformation shielded the side of target RNA (Fig. 4J and K) . In the CD conformation, the approaching of the PIWI and PAZ domains yielded the direct interactions between them, wrapping the target RNA more tightly and protecting it from solvent water (Fig. 4L and M) . To quantify this change, we calculated the solvent accessible surface area (SASA) of the target nucleotide surrounded by the PAZ, PIWI and N domains which includes the A9-C12 of target RNA in the T2 systems and its counterpart U6-A9 of target RNA in the T3 systems. The SASA value decreased from the CR state to the CL state in the T2 systems and from the CM state to the CU state in the T3 systems ( Fig. 4N and O and Supplementary Figure 10) , indicating that the closer contacts between the PAZ, PIWI and N domains induced by the miRNA extensions could protect the target RNA from water solvent. To discuss the linkage between these states to the target RNA regulation, we calculated the RMSDs between representative structures of each state and reported crystal structures (Supplementary Table 1 ). For RNA sequence and pairing pattern of target RNA and miRNA could affect the structure of AGO protein, we selected the crystal structure of seed-plus-supplementary-pairing AGO2 complex (PDB ID 6N4O) [29] as the representative structure of target inhibiting process and the crystal structure of only four central mismatches (bu4) pairing AGO2 complex (PDB ID 6NIT) [30] as the representative structure of target-directed microRNA degradation (TDMD) process to make RNA sequence and pairing pattern of the crystal structures as similar as possible to the simulation system. The RMSDs of whole protein and key domains (PAZ and PIWI domains) in AGO-RNA interaction between each representative structure and the crystal structures were calculated respectively. The result showed that RMSDs with target inhibiting conformation are smaller than that with TDMD conformation, correlated with the fact that systems we built were used to describe target inhibiting process rather than TDMD process. However, the phenomenon cannot be observed in the 23mer-T3 system. This can be ascribed to the huge difference in pairing pattern and RNA sequence between the 23mer-T3 system and the crystal structure 6N4O. The CR conformation has the lowest RMSDs with target inhibiting conformation among all states especially when calculating key domains, indicating that the crystal structure 6N4O adopted the CR conformation. This is consistent with the previous observation that CR state was the dominant conformation in the 21mer-T2 system for crystal structure 6N4O has 21mer miRNA and almost the same pairing pattern as the 21mer-T2 system. Taken together, compared to the seed-pairing-only systems, the additional supplementary pairing could induce ''closed" to ''open" motions of the AGO2 protein and promote the proximity of the PAZ and PIWI domains. Further, the proximity of the PAZ, PIWI and N domains led to the shrinkage of the supplementary chamber and tightened the mouth of the chamber, which in turn prevented the escape of target RNA from the complex and shielded it from solvent water attack. To examine how the pairing patterns and the miRNA extensions influenced the interaction of the target RNA with the AGO2 and miRNAs, we calculated the number of hydrogen bonds formed between AGO2-target RNA and miRNA-target RNA in each system. The cutoff distance and angle of hydrogen bonds are set as 3.5 Å and 135°, respectively. With the additional supplementary pairing, the number of hydrogen bonds between AGO2-target and miRNAtarget increased compared to the seed-pairing-only systems (Fig. 5A) . Additionally, the number of hydrogen bonds formed between target and protein in the T2 system was smaller than that in the T3 system while the number of hydrogen bonds formed between target and guide was just opposite. The deviation can be ascribed to the different pairing patterns and target RNA structures in the T2 and T3 systems. There are two more nucleotides pairing to miRNA at the 5 0 end of target RNA in the T3 systems (pairing to g15-g19) compared to the T2 systems (pairing to g13-g17). The additional pairing pulled the 5 0 end of target RNA from the PIWI domain of protein. Moreover, target2 has an additional 5 0 tail compared to target3 (Fig. 1D and E) , which can form additional hydrogen bonds with protein. To prove the explanation above, we calculated the hydrogen bonds formed between the PIWI domain and target RNA and between protein and the 5 0 additional tail of target RNA (A1 target , C2 target and G3 target ) in the T2 and T3 systems (Supplementary Figure 11) . We could observe that the 5 0 additional tail of target RNA in the T2 systems provided $ 6 additional hydrogen bonds with AGO2 in contrast with the T3 systems and the difference of the protein-target RNA hydrogen bonds between the T2 and T3 systems is also about 6 (Fig. 5A ). The T2 systems had about seven more hydrogen bonds between target RNA and the PIWI domain than the T3 system and the 5 0 additional tail of the T2 systems mainly formed hydrogen bonds with the PIWI domain. Based on these, we speculated that hydrogen bonds formed between the 5 0 additional tail in the T2 systems provided the most deviation in the number of target RNA-protein hydrogen bonds between the T2 and T3 systems. The number of hydrogen bonds between miRNA and target RNA in the 21mer-T3 system was slightly larger relative to the 22mer-T3 and 23mer-T3 systems. This can be explained that the release of the 3 0 end of the guide RNA relieved the conformational limitation of the PAZ domain on its tail and it could adopt a more suitable conformation to pair with the target RNA. To further elucidate the pairing between the target and guide RNAs, the pairing occupancy of each nucleotide pair was calculated using the average of occupancy of each hydrogen bond formed between the nucleotide pair. Nucleotide pairs with more than 3% variation in the pairing occupancy as the miRNA extensions were shown in Fig. 5B and C. The pairing occupancy at the edge of supplementary pairing which contains U17 miRNA -A6 target in the T2 systems and U19 miRNA -A1 target and G18 miRNA -C2 target in the T3 systems increased slightly with the miRNA extensions except for the 21mer-T3 system in which the 3 0 tail of the guide RNA released from the PAZ domain. The opening angles of the U17 miRNA -A6 target in the T2 systems and the U19 miRNA -A1 target in the T3 systems were further calculated to describe the pairing at the edge of supplementary region (Supplementary Figure 12) . With the miRNA extensions, the opening angle became smaller, suggesting a better pairing at the edge of supplementary region. The pairing at the edge of supplementary region depends on the tail of miRNA's motions towards the target RNA, which may increase the tension in the 3 0 end of miRNA. Hence, with the extended miRNA lengths, the tension can be reduced and therefore the miRNA can move together with the PAZ domain towards the target RNA. In addition to the edge pairing, the pairing occupancy of other sites seemed to vary less regularly, but generally compensated for the variation in the edge pairing caused by the miRNA extensions. This was correlated with the findings that the extended miRNA lengths did not affect the sum of hydrogen bonds between the target RNA and miRNA. Using Generalized correlation analysis, we provided an overview of the inter-residue and inter-domain correlated motions in the nine systems (see Material and Methods). The results showed that the additional supplementary pairing significantly weakened the overall correlated motions (Supplementary Figure 13A-J) . Additionally, the overall correlated movements also enhanced with the miRNA extensions. Further investigation on the interdependent motions between different domains indicated that the increased correlated motions between the PAZ domain and other domains (especially the MID and PIWI domains) were more significant in the context of overall enhanced domain connectivity with the miRNA extensions (Fig. 6A-I) . All these observations, except the 21mer-T3 system in which the 3 0 release of miRNA from the PAZ domain, reshaped the overall correlation structure and strengthened the correlated motions compared to other seed-plussupplementary-pairing systems. These results were consistent with the proximity of the PAZ, PIWI and N domains with the extended miRNA lengths, implying that miRNA extensions could enhance the domain connectivity between the PAZ, PIWI and N domains to promote the proximity of these domains to wrap around the target RNA more tightly. Community network analysis was performed to uncover the correlation network in different systems (see Material and Methods). Each community was visualized as colored circles whose area is proportional to the number of residues it contains. The strength of inter-community connection was represented by the width of sticks connecting communities. In general, communities of different systems were similar (Fig. 7A -I) and the detailed configurations of the communities in each system were shown in Supplementary Figure 14 to reveal the relationship between domains and communities. The PIWI domain consisted of Communities C, D, and F while the community F only existed in the T2 systems. The Community F diminished from the 21mer-T2 system to the 22mer-T2 system and even vanished in the 23mer-T2 system, suggesting a better motion integrity of the PIWI domain with the miRNA extensions. As for the nucleic acids, they were segmented into communities in pairs, playing a role of connecting communities along the way. In terms of where each nucleotide belonging, the 5 0 non-pairing nucleotides of target RNA in the T2 systems gradually joined Community I and moved together with miRNAs during simulations with the extension of miRNA lengths from 21mer to 23mer (Supplementary Figure 12D-F) , indicating better interactions between target and guide RNAs induced by the miRNA extensions. For the connectivity between communities, we focused on the connection between Community I (main part of the PAZ domain) and Community C (main part of the PIWI domain). In the seedpairing-only systems, the connection between the Community I and Community C was made through two additional communities which were Community G/H and Community D (Fig. 7A-C) . However, with the additional supplementary pairing, the connection was through one additional community which was Community G/H, except for the 21mer-T2 system whose connection network was influenced by the additional community F (Fig. 7D-I) . This analysis suggested that the supplementary pairing can promote the information flow between the PAZ and PIWI domains. In the seed-plus-supplementary pairing systems, the connection between Community I and C was also strengthened with the miRNA extensions through the enhancement of connections along the way while this change cannot be observed in the seed-pairingonly systems. These results indicated that the miRNA extensions can produce stronger information transduction between the PAZ and PIWI domains, which was in consistence with the generalized correlation analysis revealing that the miRNA extensions can enhance the domain connectivity between the PAZ and PIWI domains. The human AGO2 protein, acting by guiding miRNA to recognize target RNA, plays a critical role in gene expression regulation. Here, using MD simulations, we investigated how the pairing pattern and miRNA lengths influence target affinity at the atomic level, thereby providing further guidance for the design of miRNA for experimental and therapeutic use. We observed the release of the 3 0 end of the miRNA from the PAZ domain in the 21mer-T3 system, indicating that the length of the 21mer miRNA was insufficient to tolerate the increased tension created by the 3 0 tail pairing. However, when building the 21mer-T3 system, we pre-paired the target with miRNA while in real structure the target may not perfectly pair with miRNA due to excessive tension and cannot pull the miRNA from the PAZ domain. Validation of the release of 3 0 tail by experiment will be the focus of our follow-up research. Distinct pairing patterns can cause different conformational dynamics of AGO2-RNA complex. By comparing the conformational dynamics between the seed-pairing-only and seed-plussupplementary pairing systems, we confirmed the structural model that the PAZ domain shifted away from the central cleft for the opening of supplementary chamber to form additional supplementary pairing between the target RNA and miRNA. Moreover, results indicated that in addition to the supplementary guidetarget interaction, supplementary pairing could increase target affinity by promoting the proximity of the PAZ and PIWI domains to wrap target RNA more tightly. We investigated how miRNA lengths influence target affinity. Previous studies speculated that miRNA extensions increased target affinity by releasing tension in miRNA 3 0 tail to promote the supplementary interaction between guide and target. Consistent with this speculation, we found that miRNA extensions could promote pairing at the edge of supplementary region. However, the sum of hydrogen bonds formed between target and guide was similar with different miRNA lengths which can be ascribed to the compensation of inner stronger pairing. This result explained and accorded with the fact that mismatches at the edge of supplementary pairing had a minor effect on target binding. Based on this, we assumed that the enhancement of edge pairing of supplementary region was not the key factor causing the difference in target affinity with miRNA extensions. Furthermore, we put forward that miRNA extensions could induce conformational transition from ''departing" to ''approaching" with the proximity of the PAZ, PIWI and N domains. This conformational transition could tighten the mouth of the supplementary chamber, which in turn prevented the escape of target RNA from the complex and shielded it from solvent water attack. This conformational change created by the miRNA extensions can be ascribed to the enhanced domain connectivity between the PAZ, PIWI and N domains and the increased plasticity of the PAZ domain. In view of the wide experimental application of miRNA and the great potential of miRNA therapy in the treatment of tumor, cardiovascular disease, endocrine and metabolic disease, the design of miRNA mimics and anti-miRNAs has become the focus of research [11, 81, 82] . Much guidance on miRNA mimics and anti-miRNAs design has been established [83, 84] . However, the structural basis for the assembly of AGO-RNA complex needs to be further investigated for the development of new design strategies. Our study illuminated the structural basis of pairing pattern and miRNA lengths' effect on target RNA recognition, providing insights on control of pairing patterns and miRNA lengths on the design of miRNA mimics and anti-miRNAs. For the potential engineering of AGO proteins as gene silencing and even editing tools [85] , our study may provide a useful guidance to tune AGO proteins for new gene regulation tools. In summary, the collective sampling of 90 ls MD simulations revealed the effect of pairing patterns and miRNA lengths on the conformational dynamics and the assembly of AGO-RNA complexes. Compared to the seed-paring-only systems, the supple-mentary pairing can strengthen target-guide interactions and promote the proximity of the PAZ and PIWI domains, which can tighten the mouth of supplementary chamber to protect the target RNA from dissociation. Further, in the seed-plus-supplementary pairing systems, the extension of miRNA lengths could induce the increased plasticity of the PAZ domain and enhance domain connectivity between the PAZ, PIWI and N domains, and thereby promote the proximity of these domains to tighten the mouth of supplementary chamber, which in turn protected the target RNA from water solvent attacking and hindered target RNA's escaping from the AGO-RNA complexes through the mouth of supplementary chamber. Moreover, our results showed that miRNA extensions could strengthen the target-guide pairing at the edge of supplementary region while the sum of hydrogen bonds formed between target and guide RNAs were similar with different miRNA lengths and the inner stronger matching can compensate for edge weaker matching. The proximity of the PAZ, PIWI and N domains rather than the enhanced edge pairing is the key factor for the increased target affinity with the miRNA extensions. These results shed light on the structural details of the assembly of AGO-RNA complex with the supplementary pairing, which were instructive to the design of miRNA mimics and anti-miRNA for clinical and experimental use. Ethics approval and consent to participate. Not applicable. All the authors have approved and agreed to publish this manuscript. All data generated and analysed to support the conclusion of this study is included in this published manuscript and its additional supporting files. The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. The Argonaute protein family Argonaute proteins: functional insights and emerging roles Argonaute proteins: structural features, functions and emerging roles An overview of microRNAs Regulation of microRNA function in animals Analysis of the expression, function, prognosis and co-expression genes of DDX20 in gastric cancer Intestinal microbiota alterations by dietary exposure to chemicals from food cooking and processing. Application of data science for risk prediction The functions of animal microRNAs Duplexes of 21-nucleotide RNAs mediate RNA interference in cultured mammalian cells An RNA-directed nuclease mediates post-transcriptional gene silencing in Drosophila cells MicroRNA therapeutics: towards a new era for the management of cancer and other diseases miRNA nanotherapeutics for cancer Computational methods for detecting cancer hotspots Argonaute2 is the catalytic engine of mammalian RNAi The crystal structure of human argonaute2 The structure of human argonaute-2 in complex with miR-20a Argonaute divides Its RNA guide into domains with distinct functions and RNA-binding properties Structural basis for microRNA targeting Principles of microRNA-target recognition Unambiguous identification of miRNA: target site interactions by different types of ligation reactions Most mammalian mRNAs are conserved targets of microRNAs MicroRNA targeting specificity in mammals: determinants beyond seed pairing Pairing beyond the seed supports MicroRNA targeting specificity Mapping the human miRNA interactome by CLASH reveals frequent noncanonical binding An interplay of miRNA abundance and target site architecture determines miRNA activity and specificity Robust differential microRNA targeting driven by supplementary interactions in vitro Single-molecule imaging reveals that argonaute reshapes the binding properties of its nucleic acid guides Beyond the seed: structural basis for supplementary microRNA targeting Structural basis for target-directed MicroRNA degradation Non-canonical targets destabilize microRNAs in human Argonautes The ZSWIM8 ubiquitin ligase mediates target-directed microRNA degradation Target RNAdirected trimming and tailing of small silencing RNAs Masking the 5 0 terminal nucleotides of the hepatitis C virus genome by an unconventional microRNA-target RNA complex Posttranscriptional generation of miRNA variants by multiple nucleotidyl transferases contributes to miRNA transcriptome complexity Dynamic isomiR regulation in Drosophila development Computational evolution of an RNA-binding protein towards enhanced oxidized-RNA binding A systematic strategy for the investigation of vaccines and drugs targeting bacteria Programmed À1 ribosomal frameshifting from the perspective of the conformational dynamics of mRNA and ribosomes Structural insights into the mechanism of RNA recognition by the N-terminal RNAbinding domain of the SARS-CoV-2 nucleocapsid phosphoprotein Structural insight into the recognition of S-adenosyl-L-homocysteine and sinefungin in SARS Nsp16/Nsp10 RNA cap 2 0 -O-Methyltransferase Dissecting nucleotide selectivity in viral RNA polymerases Differential impact of BTK active site inhibitors on the conformational state of full-length BTK Harnessing reversed allosteric communication: a novel strategy for allosteric drug discovery Targeting a cryptic allosteric site of SIRT6 with small-molecule inhibitors that inhibit the migration of pancreatic cancer cells Mechanism of allosteric activation of SIRT6 revealed by the action of rationally designed activators Activation pathway of a G protein-coupled receptor uncovers conformational intermediates as targets for allosteric drug design Insight into the mechanism of allosteric activation of PI3Ka by oncoprotein K-Ras4B Identification of an allosteric hotspot for additive activation of PPARc in antidiabetic effects Mechanistic insights into the effect of phosphorylation on Ras conformational dynamics and its interactions with cell signaling proteins Atomic-scale insights into allosteric inhibition and evolutional rescue mechanism of Streptococcus thermophilus Cas9 by the anti-CRISPR protein AcrIIA6 Delineating the activation mechanism and conformational landscape of a class B G protein-coupled receptor glucagon receptor Untangling dual-targeting therapeutic mechanism of epidermal growth factor receptor (EGFR) based on reversed allosteric communication Markov state models and molecular dynamics simulations provide understanding of the nucleotide-dependent dimerization-based activation of lrrk2 roc domain Allosteric methods and their applications: facilitating the discovery of allosteric drugs and the investigation of allosteric mechanisms Small molecule allosteric modulators of G-protein-coupled receptors: drug-target interactions Allosteric modulator discovery: from serendipity to structure-based design Insights into the effect of molecular crowding on the structure, interactions and functions of siRNA-PAZ complex through molecular dynamics studies Markov state models reveal a two-step mechanism of miRNA loading into the human argonaute protein: selective binding followed by structural re-arrangement Probing the binding interactions between chemically modified siRNAs and human argonaute 2 using microsecond molecular dynamics simulations Two symmetric arginine residues play distinct roles in thermus thermophilus argonaute DNA guide strand-mediated DNA target cleavage Prediction of the miRNA interactome -established methods and upcoming perspectives Prediction methods for microRNA targets in bilaterian animals: toward a better understanding by biologists Watermediated recognition of t1-adenosine anchors Argonaute2 to microRNA targets Refinement of the Cornell et al. Nucleic acids force field based on reference quantum chemical calculations of glycosidic torsion profiles Refinement of the AMBER force field for nucleic acids: improving the description of a/c conformers ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB Comparison of simple potential functions for simulating liquid water Atomic-level accuracy in simulations of large protein crystals Integrated analysis of copy number variation and genome-wide expression profiling in colorectal cancer tissues Essential dynamics: foundation and applications Generalized correlation for biomolecular dynamics 0: a package for molecular simulation and trajectory analysis Network view: 3D display and analysis of protein ÁrNA interaction networks Community structure in social and biological networks Markov state models: from an art to a science Computational methods for protein localization prediction A review of computational tools for generating metagenome-assembled genomes from metagenomic sequencing data Discovery of cryptic allosteric sites using reversed allosteric communication by a combined computational and experimental strategy Deactivation pathway of Ras GTPase underlies conformational substates as targets for drug design MicroRNA-based therapeutics in cardiovascular disease: screening and delivery to the target Developing MicroRNA therapeutics Anti-miRNA oligonucleotides: a comprehensive guide for design Guidelines for the optimal design of miRNA-based shRNAs Strategies for mitochondrial gene editing This work was supported by the National Natural Science Foundation of China (82172846 and 81901255). Supplementary data to this article can be found online at https://doi.org/10.1016/j.csbj.2022.03.001.