We present the color and flux variability analysis at 3.4 µm (W 1-band) and 4.6 µm (W 2-band) of 492 narrow-line Seyfert 1 (NLSy1) galaxies using archival data from the Wide-field Infrared Survey Explorer (WISE). In the WISE color-color, (W 1 − W 2) versus (W 2 − W 3) diagram, ∼ 58% of the NLSy1 galaxies of our sample lie in the region occupied by the blazar category of active galactic nuclei (AGN). The mean W 1 − W 2 color of candidate variable NLSy1 galaxies is 0.99 ± 0.18 mag. The average amplitude of variability is 0.11 ± 0.07 mag in long-term (multi-year) with no difference in variability between W 1 and W 2-bands. The W 1 − W 2 color of NLSy1 galaxies is anti-correlated with the relative strength of [O III] to Hβ, strongly correlated with continuum luminosity, black hole mass, and Eddington ratio. The long-term amplitude of variability shows weak anti-correlation with the Fe II strength, continuum luminosity and Eddington ratio. A positive correlation between color as well as the amplitude of variability with the radio power at 1.4 GHz was found for the radio-detected NLSy1 galaxies. This suggests non-thermal synchrotron contribution to the mid-infrared color and flux variability in radio-detected NLSy1 galaxies.