This study focuses on the analysis of mid-infrared variability in a sample of high-mass young stellar objects (YSOs) associated with the cataloged sources from the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL). The Near-Earth Object Wide-field Infrared Survey Explorer Reactivation Mission (NEOWISE) database was used to explore the long-term mid-infrared variability of these high-mass YSOs at a half-year scale. After matching with NEOWISE photometric measurements, a total of 2230 ATLASGAL sources were selected for the variability analysis, out of which 717 were identified as variables. The derived proportions of variables at different evolutionary stages show that the variability rate of high-mass YSOs is highest during the YSO stage and decreases with evolution toward the H ii region stage, resembling the behavior of low-mass YSOs. The variables can be classified into six types based on their light curves, divided into two categories: secular (linear, sin, sin+linear) and stochastic variables (burst, drop, and irregular). The magnitude-color variations observed in ~160 secular variables can be mainly divided into "bluer when brighter/redder when dimming" and "redder when brighter/bluer when dimming", likely originating from changes in accretion rate or the effect of extinction due to obscuration. Moreover, several episodic accretion candidates were selected for further observational studies.