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.