The low-lying excitation observables of even-even nuclei, e.g., the excitation energies of the first $2^{+}$ and $4^{+}$ states, $E(2_1^{+})$ and $E(4_1^{+})$, and the reduced electric quadrupole transition probability between the ground state and the first $2^{+}$ state $B(E2;0_1^{+}\to 2_1^{+})$, serve as promising experimental probes for nuclear shell evolution in unstable nuclei far from $\beta$-stability line. Previous machine-learning studies on low-lying nuclear states have typically employed independent single-task models that neglect the intrinsic correlations among $B(E2;0_1^{+}\to 2_1^{+})$, $E(2_1^{+})$, and $E(4_1^{+})$, potentially leading to unphysical divergences in practical extrapolation. In addition, the manual selection of hyperparameters introduces subjective tuning bias and model dependence. To overcome these limitations, we propose a multi-task Gaussian process (MTGP) approach that jointly predicts the three observables while preserving their mutual physical couplings. The MTGP hyperparameters are automatically optimized via Bayesian optimization. The examination of the established MTGP method demonstrates its significant improvement in predictive robustness over previous single-task approaches, both in terms of the prediction accuracy and the description of inherent correlations among the three observables. We then apply the MTGP model to predict unmeasured $B(E2;0_1^{+}\to 2_1^{+})$, $E(2_1^{+})$, and $E(4_1^{+})$ values around the magic shells $N,Z=20$ and $28$. Systematic analysis of the MTGP-predicted observables along their respective isotonic and isotopic chains reveal nuanced shell-evolution patterns at $N,Z=20$ and $28$: the $N=20$ closure weakens in the neutron-rich $Z=8-12$ region, strengthens at $Z=14-20$, and weakens again for $Z\ge 22$; the $N=28$ closure is quenched at $Z=12,14$, robust at $Z=16-28$, and shows signs of weakening in ${}^{58}$Zn and ${}^{60}$Ge; the $Z=20$ closure persists in ${}^{34,36}$Ca but gradually weakens beyond $N=32$; the $Z=28$ closure remains strong in ${}^{48,50,52}$Ni and weakens in neutron-rich nickel isotopes beyond $N=46$. These data-driven results support that the MTGP-based systematics capture the evolution of shell gaps and collectivity, suggesting MTGP as a potential tool for nuclear structure studies.