We present the first lattice QCD calculation of the next-to-leading order (NLO) hadronic vacuum polarization (HVP) contribution to the muon anomalous magnetic moment with sub-percent precision. We employ the time-momentum representation combined with the spatially summed vector correlator computed on CLS ensembles with $N_{\mathrm{f}}=2+1$ flavors of $\mathrm{O}(a)$-improved Wilson fermions, spanning six lattice spacings ($0.039$-$0.097\,$fm) and a range of pion masses including the physical value. After accounting for finite-size corrections and isospin-breaking effects, we obtain in the continuum limit $a_\mu^{\mathrm{hvp,\,nlo}} = (-101.57 \pm 0.26_{\rm stat} \pm 0.54_{\rm syst}) \times 10^{-11}$, corresponding to a total relative error of $0.6\%$. Our result lies $1.4\sigma$ below the estimate of the 2025 White Paper update and is two times more precise. It also shows a tension of $4.6\sigma$ with data-driven evaluations based on hadronic cross section measurements prior to the CMD-3 result.