Modeling loop reorganization free energies of acetylcholinesterase: a comparison of explicit and implicit solvent models

The treatment of hydration effects in protein dynamics simulations varies in model complexity and spans the range from the computationally intensive microscopic evaluation to simple dielectric screening of charge-charge interactions. This paper compares different solvent models applied to the problem of estimating the free-energy difference between two loop conformations in acetylcholinesterase. Molecular dynamics (MD) simulations were used to sample potential energy surfaces of the two basins with solvent treated by means of explicit and implicit methods. Implicit solvent methods studied include the generalized Born (GB) model, atomic solvation potential (ASP), and the distance-dependent dieletric constant. By using the linear response approximation (LRA), the explicit solvent calculations determined a free-energy difference that is in excellent agreement with the experimental estimate, while rescoring the protein conformations with GB or the Poisson equation showed inconsistent and inferior results. While the approach of rescoring conformations from explicit water simulations with implicit solvent models is popular among many applications, it perturbs the energy landscape by changing the solvent contribution to microstates without conformational relaxation, thus leading to non-optimal solvation free energies. Calculations applying MD with a GB solvent model produced results of comparable accuracy as observed with LRA, yet the electrostatic free-energy terms were significantly different due to optimization on a potential energy surface favored by an implicit solvent reaction field. The simpler methods of ASP and the distance-dependent scaling of the dielectric constant both produced considerable distortions in the protein internal free-energy terms and are consequently unreliable.     

Hierarchical Cascades of Instability Govern the Mechanics of Coiled Coils: Helix Unfolding Precedes Coil Unzipping

Coiled coils are a fundamental emergent motif in proteins found in structural biomaterials, consisting of α-helical secondary structures wrapped in a supercoil. A fundamental question regarding the thermal and mechanical stability of coiled coils in extreme environments is the sequence of events leading to the disassembly of individual oligomers from the universal coiled-coil motifs. To shed light on this phenomenon, here we report atomistic simulations of a trimeric coiled coil in an explicit water solvent and investigate the mechanisms underlying helix unfolding and coil unzipping in the assembly. We employ advanced sampling techniques involving steered molecular dynamics and metadynamics simulations to obtain the free-energy landscapes of single-strand unfolding and unzipping in a three-stranded assembly. Our comparative analysis of the free-energy landscapes of instability pathways shows that coil unzipping is a sequential process involving multiple intermediates. At each intermediate state, one heptad repeat of the coiled coil first unfolds and then unzips due to the loss of contacts with the hydrophobic core. This observation suggests that helix unfolding facilitates the initiation of coiled-coil disassembly, which is confirmed by our 2D metadynamics simulations showing that unzipping of one strand requires less energy in the unfolded state compared with the folded state. Our results explain recent experimental findings and lay the groundwork for studying the hierarchical molecular mechanisms that underpin the thermomechanical stability/instability of coiled coils and similar protein assemblies.     

The investigation of the secondary structure propensities and free-energy landscapes of peptide ligands by replica exchange molecular dynamics simulations

The conformational states of two peptide sequences that bind to staphylococcal enterotoxin B are sampled by replica exchange molecular dynamic (REMD) simulations in explicit water. REMD simulations were treated with 52 replicas in the range of 280–501 K for both peptides. The conformational ensembles of both peptides are dominated by random coil, bend and turn structures with a small amount of helical structures for each temperature. In addition, while an insignificant presence of β-bridge structures were observed for both peptides, the β-sheet structure was observed only for peptide 3. The results obtained from simulations at 300 K are consistent with the experimental results obtained from circular dichroism spectroscopy. From the analysis of REMD results, we also calculated hydrophobic and hydrophilic solvent accessible surface areas for both peptides, and it was observed that the hydrophobic segments of the peptides tend to form bend or turn structures. Moreover, the free-energy landscapes of both peptides were obtained by principal component analysis to understand how the secondary structural properties change according to their complex space. From the free-energy analysis, we have found several minima for both peptides at decreased temperature. For these obvious minima of both peptides, it was observed that the random coil, bend and turn structures are still dominant and the helix, β-bridge or β-sheet structures can appear or disappear with respect to minima. On the other hand, when we compare the results of REMD with conventional MD simulations for these peptides, the configurations of peptide 3 might be trapped in energy minima during the conventional MD simulations. Hence, it can be said that the REMD simulations have provided a sufficiently high sampling efficiency.     

The Free Energy Profile of Tubulin Straight-Bent Conformational Changes,with Implications for Microtubule Assembly and Drug Discovery

αβ-tubulin dimers need to convert between a ‘bent’ conformation observed for free dimers in solution and a ‘straight’ conformation required for incorporation into the microtubule lattice. Here, we investigate the free energy landscape of αβ-tubulin using molecular dynamics simulations, emphasizing implications for models of assembly, and modulation of the conformational landscape by colchicine, a tubulin-binding drug that inhibits microtubule polymerization. Specifically, we performed molecular dynamics, potential-of-mean force simulations to obtain the free energy profile for unpolymerized GDP-bound tubulin as a function of the ∼12° intradimer rotation differentiating the straight and bent conformers. Our results predict that the unassembled GDP-tubulin heterodimer exists in a continuum of conformations ranging between straight and bent, but, in agreement with existing structural data, suggests that an intermediate bent state has a lower free energy (by ∼1 kcal/mol) and thus dominates in solution. In agreement with predictions of the lattice model of microtubule assembly, lateral binding of two αβ-tubulins strongly shifts the conformational equilibrium towards the straight state, which is then ∼1 kcal/mol lower in free energy than the bent state. Finally, calculations of colchicine binding to a single αβ-tubulin dimer strongly shifts the equilibrium toward the bent states, and disfavors the straight state to the extent that it is no longer thermodynamically populated.     

Ab initio folding simulation of Trpcage by replica exchange with hybrid Hamiltonian

Replica-exchange molecular dynamics simulations with hybrid Hamiltonian in explicit solvent were performed to study the folding of a designed 20-residue miniprotein, Trpcage, from a fully extended structure. During the simulations several folding/unfolding events happened. In the folded states the majority of experimentally observed NMR NOE restraints are satisfied. The folded structures have root mean squared deviation of 2.0 A with respect to the NMR structures considering all heavy atoms. The free-energy surface constructed by the hybrid Hamiltonian simulations is similar to the one built by a standard replica-exchange simulation which started from the native structure. Consistent with previous experimental observation, a pre-existing hydrophobic collapse in the unfolded state is detected by investigating the desolvation behavior of Trpcage. At room temperature, an intermediate state featured by a misfolded core, a nearly formed alpha-helix segment and an absence of 3(10)-helix is found. The replica exchange with hybrid Hamiltonian method is shown here to be capable of resolving the folding picture of the miniprotein.     

K+ Block Is the Mechanism of Functional Asymmetry in Bacterial Nav Channels

Crystal structures of several bacterial Na_v channels have been recently published and molecular dynamics simulations of ion permeation through these channels are consistent with many electrophysiological properties of eukaryotic channels. Bacterial Na_v channels have been characterized as functionally asymmetric, and the mechanism of this asymmetry has not been clearly understood. To address this question, we combined non-equilibrium simulation data with two-dimensional equilibrium unperturbed landscapes generated by umbrella sampling and Weighted Histogram Analysis Methods for multiple ions traversing the selectivity filter of bacterial Na_vAb channel. This approach provided new insight into the mechanism of selective ion permeation in bacterial Na_v channels. The non-equilibrium simulations indicate that two or three extracellular K⁺ ions can block the entrance to the selectivity filter of Na_vAb in the presence of applied forces in the inward direction, but not in the outward direction. The block state occurs in an unstable local minimum of the equilibrium unperturbed free-energy landscape of two K⁺ ions that can be ‘locked’ in place by modest applied forces. In contrast to K⁺, three Na⁺ ions move favorably through the selectivity filter together as a unit in a loose “knock-on” mechanism of permeation in both inward and outward directions, and there is no similar local minimum in the two-dimensional free-energy landscape of two Na⁺ ions for a block state. The useful work predicted by the non-equilibrium simulations that is required to break the K⁺ block is equivalent to large applied potentials experimentally measured for two bacterial Na_v channels to induce inward currents of K⁺ ions. These results illustrate how inclusion of non-equilibrium factors in the simulations can provide detailed information about mechanisms of ion selectivity that is missing from mechanisms derived from either crystal structures or equilibrium unperturbed free-energy landscapes.     

Basis for the Specificity and Activation of the Serpin Protein Z-dependent Proteinase Inhibitor (ZPI) as an Inhibitor of Membrane-associated Factor Xa

The serpin ZPI is a protein Z (PZ)-dependent specific inhibitor of membrane-associated factor Xa (fXa) despite having an unfavorable P1 Tyr. PZ accelerates the inhibition reaction ∼2000-fold in the presence of phospholipid and Ca²⁺. To elucidate the role of PZ, we determined the x-ray structure of Gla-domainless PZ (PZ_ΔGD) complexed with protein Z-dependent proteinase inhibitor (ZPI). The PZ pseudocatalytic domain bound ZPI at a novel site through ionic and polar interactions. Mutation of four ZPI contact residues eliminated PZ binding and membrane-dependent PZ acceleration of fXa inhibition. Modeling of the ternary Michaelis complex implicated ZPI residues Glu-313 and Glu-383 in fXa binding. Mutagenesis established that only Glu-313 is important, contributing ∼5–10-fold to rate acceleration of fXa and fXIa inhibition. Limited conformational change in ZPI resulted from PZ binding, which contributed only ∼2-fold to rate enhancement. Instead, template bridging from membrane association, together with previously demonstrated interaction of the fXa and ZPI Gla domains, resulted in an additional ∼1000-fold rate enhancement. To understand why ZPI has P1 tyrosine, we examined a P1 Arg variant. This reacted at a diffusion-limited rate with fXa, even without PZ, and predominantly as substrate, reflecting both rapid acylation and deacylation. P1 tyrosine thus ensures that reaction with fXa or most other arginine-specific proteinases is insignificant unless PZ binds and localizes ZPI and fXa on the membrane, where the combined effects of Gla-Gla interaction, template bridging, and interaction of fXa with Glu-313 overcome the unfavorability of P1 Tyr and ensure a high rate of reaction as an inhibitor.     

Small Peptide Binding Stiffens the Ubiquitin-like Protein SUMO1

Posttranslational modification by small ubiquitin-like modifiers (SUMOs), known as SUMOylation, is a key regulatory event in many eukaryotic cellular processes in which SUMOs interact with a large number of target proteins. SUMO binding motifs (SBMs) are small peptides derived from these target proteins that interact noncovalently with SUMOs and induce conformational changes. To determine the effect of SBMs on the mechanical properties of SUMO1 (the first member of the human SUMO family), we performed single-molecule force spectroscopy experiments on SUMO1/SBM complexes. The unfolding force of SUMO1 (at a pulling speed of 400 nm/s) increased from ∼130 pN to ∼170 pN upon binding to SBMs, indicating mechanical stabilization upon complexation. Pulling-speed-dependent experiments and Monte Carlo simulations measured a large decrease in distance to the unfolding transition state for SUMO1 upon SBM binding, which is by far the largest change measured for any ligand binding protein. The stiffness of SUMO1 (measured as a spring constant for the deformation response along the line joining the N- and C-termini) increased upon SBM binding from ∼1 N/m to ∼3.5 N/m. The relatively higher flexibility of ligand-free SUMO1 might play a role in accessing various conformations before binding to a target.     

Short divalent ethacrynic amides as pro-inhibitors of glutathione S-transferase isozyme Mu and potent sensitisers of cisplatin-resistant ovarian cancers

The linking of ethacrynic acid with ethylenediamine and 1,4-butanediamine gave EDEA and BDEA, respectively, as membrane-permeable divalent pro-inhibitors of glutathione S-transferase (GST). Their divalent glutathione conjugates showed subnanomolar inhibition and divalence-binding to GSTmu (GSTM) (PDB: 5HWL) at ∼0.35 min⁻¹. In cisplatin-resistant SK-OV-3, COC1, SGC7901 and A549 cells, GSTM activities probed by 15 nM BDEA or EDEA revealed 5-fold and 1.0-fold increases in cisplatin-resistant SK-OV-3 and COC1 cells, respectively, in comparison with the susceptible parental cells. Being tolerable by HEK293 and LO2 cells, BDEA at 0.2 μM sensitised resistant SK-OV-3 and COC1 cells by ∼3- and ∼5-folds, respectively, released cytochrome c and increased apoptosis; EDEA at 1.0 μM sensitised resistant SK-OV-3 and A549 cells by ∼5- and ∼7-fold, respectively. EDEA at 1.7 μg/g sensitised resistant SK-OV-3 cells to cisplatin at 3.3 μg/g in nude mouse xenograft model. BDEA and EDEA are promising leads for probing cellular GSTM and sensitising cisplatin-resistant ovarian cancers.     

Weight average approaches for predicting dynamical properties of biomolecules

Recent advances in atomistic molecular dynamics (MD) simulations of biomolecules allow us to explore their conformational spaces widely, observing large-scale conformational fluctuations or transitions between distinct structures. To reproduce or refine experimental data using MD simulations, structure ensembles, which are characterized by multiple structures and their statistical weights on the rugged free-energy landscapes, are often used. Here, we summarize weight average approaches for various experimental measurements. Weight average approaches are now applied to hybrid quantum mechanics/molecular mechanics MD simulations to predict fast vibrational motions in a protein with a high accuracy for better understanding of molecular functions from atomic structures.     

Spatially resolved free-energy contributions of native fold and molten-globule-like Crambin

The folding stability of a protein is governed by the free-energy difference between its folded and unfolded states, which results from a delicate balance of much larger but almost compensating enthalpic and entropic contributions. The balance can therefore easily be shifted by an external disturbance, such as a mutation of a single amino acid or a change of temperature, in which case the protein unfolds. Effects such as cold denaturation, in which a protein unfolds because of cooling, provide evidence that proteins are strongly stabilized by the solvent entropy contribution to the free-energy balance. However, the molecular mechanisms behind this solvent-driven stability, their quantitative contribution in relation to other free-energy contributions, and how the involved solvent thermodynamics is affected by individual amino acids are largely unclear. Therefore, we addressed these questions using atomistic molecular dynamics simulations of the small protein Crambin in its native fold and a molten-globule-like conformation, which here served as a model for the unfolded state. The free-energy difference between these conformations was decomposed into enthalpic and entropic contributions from the protein and spatially resolved solvent contributions using the nonparametric method Per|Mut. From the spatial resolution, we quantified the local effects on the solvent free-energy difference at each amino acid and identified dependencies of the local enthalpy and entropy on the protein curvature. We identified a strong stabilization of the native fold by almost 500 kJ mol⁻¹ due to the solvent entropy, revealing it as an essential contribution to the total free-energy difference of (53 ± 84) kJ mol⁻¹. Remarkably, more than half of the solvent entropy contribution arose from induced water correlations.     

Enhanced sampling simulations to construct free-energy landscape of protein–partner substrate interaction

Molecular dynamics (MD) simulations using all-atom and explicit solvent models provide valuable information on the detailed behavior of protein–partner substrate binding at the atomic level. As the power of computational resources increase, MD simulations are being used more widely and easily. However, it is still difficult to investigate the thermodynamic properties of protein–partner substrate binding and protein folding with conventional MD simulations. Enhanced sampling methods have been developed to sample conformations that reflect equilibrium conditions in a more efficient manner than conventional MD simulations, thereby allowing the construction of accurate free-energy landscapes. In this review, we discuss these enhanced sampling methods using a series of case-by-case examples. In particular, we review enhanced sampling methods conforming to trivial trajectory parallelization, virtual-system coupled multicanonical MD, and adaptive lambda square dynamics. These methods have been recently developed based on the existing method of multicanonical MD simulation. Their applications are reviewed with an emphasis on describing their practical implementation. In our concluding remarks we explore extensions of the enhanced sampling methods that may allow for even more efficient sampling.     

The Allosteric Mechanism of Activation of Antithrombin as an Inhibitor of Factor IXa and Factor Xa: HEPARIN-INDEPENDENT FULL ACTIVATION THROUGH MUTATIONS ADJACENT TO HELIX D*

Allosteric conformational changes in antithrombin induced by binding a specific heparin pentasaccharide result in very large increases in the rates of inhibition of factors IXa and Xa but not of thrombin. These are accompanied by CD, fluorescence, and NMR spectroscopic changes. X-ray structures show that heparin binding results in extension of helix D in the region 131–136 with coincident, and possibly coupled, expulsion of the hinge of the reactive center loop. To examine the importance of helix D extension, we have introduced strong helix-promoting mutations in the 131–136 region of antithrombin (YRKAQK to LEEAAE). The resulting variant has endogenous fluorescence indistinguishable from WT antithrombin yet, in the absence of heparin, shows massive enhancements in rates of inhibition of factors IXa and Xa (114- and 110-fold, respectively), but not of thrombin, together with changes in near- and far-UV CD and ¹H NMR spectra. Heparin binding gives only ∼3–4-fold further rate enhancement but increases tryptophan fluorescence by ∼23% without major additional CD or NMR changes. Variants with subsets of these mutations show intermediate activation in the absence of heparin, again with basal fluorescence similar to WT and large increases upon heparin binding. These findings suggest that in WT antithrombin there are two major complementary sources of conformational activation of antithrombin, probably involving altered contacts of side chains of Tyr-131 and Ala-134 with core hydrophobic residues, whereas the reactive center loop hinge expulsion plays only a minor additional role.     

Conformational energy landscape of the acyl pocket loop in acetylcholinesterase: a Monte Carlo-generalized Born model study

The X-ray crystal structure of the reaction product of acetylcholinesterase (AChE) with the inhibitor diisopropylphosphorofluoridate (DFP) showed significant structural displacement in a loop segment of residues 287-290. To understand this conformational selection, a Monte Carlo (MC) simulation study was performed of the energy landscape for the loop segment. A computational strategy was applied by using a combined simulated annealing and room temperature Metropolis sampling approach with solvent polarization modeled by a generalized Born (GB) approximation. Results from thermal annealing reveal a landscape topology of broader basin opening and greater distribution of energies for the displaced loop conformation, while the ensemble average of conformations at 298 K favored a shift in populations toward the native by a free-energy difference in good agreement with the estimated experimental value. Residue motions along a reaction profile of loop conformational reorganization are proposed where Arg-289 is critical in determining electrostatic effects of solvent interaction versus Coulombic charging.     

A Conformational Sampling Model for Radical Catalysis in Pyridoxal Phosphate- and Cobalamin-dependent Enzymes

Cobalamin-dependent enzymes enhance the rate of C–Co bond cleavage by up to ∼10¹²-fold to generate cob(II)alamin and a transient adenosyl radical. In the case of the pyridoxal 5′-phosphate (PLP) and cobalamin-dependent enzymes lysine 5,6-aminomutase and ornithine 4,5 aminomutase (OAM), it has been proposed that a large scale domain reorientation of the cobalamin-binding domain is linked to radical catalysis. Here, OAM variants were designed to perturb the interface between the cobalamin-binding domain and the PLP-binding TIM barrel domain. Steady-state and single turnover kinetic studies of these variants, combined with pulsed electron-electron double resonance measurements of spin-labeled OAM were used to provide direct evidence for a dynamic interface between the cobalamin and PLP-binding domains. Our data suggest that following ligand binding-induced cleavage of the Lys⁶²⁹-PLP covalent bond, dynamic motion of the cobalamin-binding domain leads to conformational sampling of the available space. This supports radical catalysis through transient formation of a catalytically competent active state. Crucially, it appears that the formation of the state containing both a substrate/product radical and Co(II) does not restrict cobalamin domain motion. A similar conformational sampling mechanism has been proposed to support rapid electron transfer in a number of dynamic redox systems.     

Asialoglycoprotein receptor 1 is a novel PCSK9-independent ligand of liver LDLR cleaved by furin

The hepatic carbohydrate-recognizing asialoglycoprotein receptor (ASGR1) mediates the endocytosis/lysosomal degradation of desialylated glycoproteins following binding to terminal galactose/N-acetylgalactosamine. Human heterozygote carriers of ASGR1 deletions exhibit ∼34% lower risk of coronary artery disease and ∼10% to 14% reduction of non-HDL cholesterol. Since the proprotein convertase PCSK9 is a major degrader of the low-density lipoprotein receptor (LDLR), we investigated the degradation and functionality of LDLR and/or PCSK9 by endogenous/overexpressed ASGR1 using Western blot and immunofluorescence in HepG2-naïve and HepG2-PCSK9-knockout cells. ASGR1, like PCSK9, targets LDLR, and both independently interact with/enhance the degradation of the receptor. This lack of cooperativity between PCSK9 and ASGR1 was confirmed in livers of wildtype (WT) and Pcsk9^−/− mice. ASGR1 knockdown in HepG2-naïve cells significantly increased total (∼1.2-fold) and cell-surface (∼4-fold) LDLR protein. In HepG2-PCSK9-knockout cells, ASGR1 silencing led to ∼2-fold higher levels of LDLR protein and DiI (1,1′-dioctadecyl-3,3,3′,3′-tetramethylindocarbocyanine perchlorate)-LDL uptake associated with ∼9-fold increased cell-surface LDLR. Overexpression of WT-ASGR1/2 primarily reduced levels of immature non-O-glycosylated LDLR (∼110 kDa), whereas the triple Ala-mutant of Gln240/Trp244/Glu253 (characterized by loss of carbohydrate binding) reduced expression of the mature form of LDLR (∼150 kDa), suggesting that ASGR1 binds the LDLR in both a sugar-dependent and -independent fashion. The protease furin cleaves ASGR1 at the RKMK¹⁰³↓ motif into a secreted form, likely resulting in a loss of function on LDLR. Altogether, we demonstrate that LDLR is the first example of a liver-receptor ligand of ASGR1. We conclude that silencing of ASGR1 and PCSK9 may lead to higher LDL uptake by hepatocytes, thereby providing a novel approach to further reduce LDL cholesterol levels.     

Hierarchical Cascades of Instability Govern the Mechanics of Coiled Coils: Helix Unfolding Precedes Coil Unzipping

Coiled coils are a fundamental emergent motif in proteins found in structural biomaterials, consisting of α-helical secondary structures wrapped in a supercoil. A fundamental question regarding the thermal and mechanical stability of coiled coils in extreme environments is the sequence of events leading to the disassembly of individual oligomers from the universal coiled-coil motifs. To shed light on this phenomenon, here we report atomistic simulations of a trimeric coiled coil in an explicit water solvent and investigate the mechanisms underlying helix unfolding and coil unzipping in the assembly. We employ advanced sampling techniques involving steered molecular dynamics and metadynamics simulations to obtain the free-energy landscapes of single-strand unfolding and unzipping in a three-stranded assembly. Our comparative analysis of the free-energy landscapes of instability pathways shows that coil unzipping is a sequential process involving multiple intermediates. At each intermediate state, one heptad repeat of the coiled coil first unfolds and then unzips due to the loss of contacts with the hydrophobic core. This observation suggests that helix unfolding facilitates the initiation of coiled-coil disassembly, which is confirmed by our 2D metadynamics simulations showing that unzipping of one strand requires less energy in the unfolded state compared with the folded state. Our results explain recent experimental findings and lay the groundwork for studying the hierarchical molecular mechanisms that underpin the thermomechanical stability/instability of coiled coils and similar protein assemblies.     

Comprehensive Characterization of Human Genome Variation by High Coverage Whole-Genome Sequencing of Forty Four Caucasians

Whole genome sequencing studies are essential to obtain a comprehensive understanding of the vast pattern of human genomic variations. Here we report the results of a high-coverage whole genome sequencing study for 44 unrelated healthy Caucasian adults, each sequenced to over 50-fold coverage (averaging 65.8×). We identified approximately 11 million single nucleotide polymorphisms (SNPs), 2.8 million short insertions and deletions, and over 500,000 block substitutions. We showed that, although previous studies, including the 1000 Genomes Project Phase 1 study, have catalogued the vast majority of common SNPs, many of the low-frequency and rare variants remain undiscovered. For instance, approximately 1.4 million SNPs and 1.3 million short indels that we found were novel to both the dbSNP and the 1000 Genomes Project Phase 1 data sets, and the majority of which (∼96%) have a minor allele frequency less than 5%. On average, each individual genome carried ∼3.3 million SNPs and ∼492,000 indels/block substitutions, including approximately 179 variants that were predicted to cause loss of function of the gene products. Moreover, each individual genome carried an average of 44 such loss-of-function variants in a homozygous state, which would completely “knock out” the corresponding genes. Across all the 44 genomes, a total of 182 genes were “knocked-out” in at least one individual genome, among which 46 genes were “knocked out” in over 30% of our samples, suggesting that a number of genes are commonly “knocked-out” in general populations. Gene ontology analysis suggested that these commonly “knocked-out” genes are enriched in biological process related to antigen processing and immune response. Our results contribute towards a comprehensive characterization of human genomic variation, especially for less-common and rare variants, and provide an invaluable resource for future genetic studies of human variation and diseases.