We report numerically accurate path integral Monte Carlo results for harmonically confined two-dimensional quantum dots containing up to N=60 interacting electrons. The finite-temperature values are extrapolated to 0 K and zero time step in order to provide precise upper-bound energies. The ground-state energies are compared against coupled-cluster and diffusion Monte Carlo results available in the literature for N≤20. We also provide Padé fits for the energies as a function of N for different strengths of the confining potential. The fits deviate less than 0.25% from the path integral Monte Carlo data. Overall, our upper-bound estimates for the ground-state energies have lower values than previous diffusion Monte Carlo benchmarks due to the accurate nodal surface in our simulations. Hence, our results set a new numerical benchmark for two-dimensional (spin-unpolarized) quantum dots up to a large number of electrons.