The JUNO first data and the recent neutrino global fit results are implemented in the sum rule from the residual ${\mathbb{Z}}_2^s$ and ${\bar{\mathbb{Z}}}_2^s$ symmetries to make prediction of the leptonic Dirac CP phase δD. Without involving model parameters, the probability distribution of δD can be readily obtained from the experimental measurements of the three mixing angles. We then confront the theoretical predictions with the global fit results for the CP phase as well as the T2K and NOvA joint analysis for their CP measurements to give the data preference of the two residual symmetries with Bayes factor for both normal and inverted orderings. We further extend our analysis to a two-dimensional probability distribution to fully explore the correlation between the CP phase δD and the atmospheric angle θa ≡  θ 23.