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 nth 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.