  We derive higher moments in the statistical distribution of OPE coefficients
in holographic 2D CFTs, and show that such moments correspond to multiboundary
Euclidean wormholes in pure 3D gravity. The n-th cyclic non-Gaussian
contraction of heavy-heavy-light OPE coefficients follows from crossing
symmetry of the thermal n-point function. We derive universal expressions for
the cubic and quartic moments and demonstrate that their scaling with the
microcanonical entropy agrees with a generalization of the Eigenstate
Thermalization Hypothesis. Motivated by this result, we conjecture that the
full statistical ensemble of OPE data is fixed by three premises: typicality,
crossing symmetry and modular invariance. Together, these properties give
predictions for non-factorizing observables, such as the generalized spectral
form factor. Using the Virasoro TQFT, we match these connected averages to new
on-shell wormhole topologies with multiple boundary components. Lastly, we
study and clarify examples where the statistics of heavy operators are not
universal and depend on the light operator spectrum. We give a gravitational
interpretation to these corrections in terms of Wilson loops winding around
non-trivial cycles in the bulk.
