A key goal in molecular evolution is to understand why one evolutionary trajectory is taken rather than others. One important determinant of these outcomes is epistasis, where the effect of mutation depends on the presence or absence of other mutations. While pairwise epistasis has been studied extensively, much less is known about high-order epistasis between three or more mutations. If present, high-order interactions could lead to a profound memory effect in evolution, with early mutations strongly shaping evolutionary outcomes. To investigate high-order epistasis, we analyzed a series of published, experimental genotype-phenotype maps using a robust statistical model. We found extensive high-order epistasis, with statistically-significant interactions between up to six mutations. Removing these interactions dramatically altered outcomes of evolutionary simulations in these maps, revealing that high-order epistasis can indeed shape evolutionary outcomes. We next investigated the origins of this epistasis, finding that the epistasis in some datasets could be explained by a globally nonlinear genotype-phenotype map. In others, it results from specific interactions between mutations. Using a biophysical model of proteins, we then showed that the observed patterns of high-order epistasis arise naturally from the ensemble nature of biomolecular systems. From these results, we propose that high-order epistasis is the rule rather than the exception in molecular evolution, and that this is a natural consequence of the physical properties of biomolecular systems. This implies that, in general, the effect of mutations will be different if they occur early or late in evolution, and that early evolutionary events strongly constrain future ones.