  We present the final results from a multi-year study of the in-medium
spectral properties of heavy quarkonium bound states on the lattice. In this
work we combine high statistics $N_f=2+1$ ensembles from the HotQCD
collaboration with the effective theory NRQCD and improved Bayesian spectral
reconstruction methods. We corroborate with high precision the hierarchical
in-medium modification of quarkonium states with respect to their vacuum
binding energy and provide updated values on melting temperatures. In
particular we are able to understand previous disagreements between different
Bayesian methods as resulting from underestimated systematic uncertainties. The
main quantitative result is a robust determination of the in-medium mass shifts
of quarkonium ground states, which we find are negative, consistent with the
behavior observed in strongly coupled pNRQCD potential based computations.
