All neutron star progenitors in neutron-star High-Mass X-ray Binaries (NS HMXBs) undergo a supernova event that may lead to a significant natal kick impacting the motion of the whole binary system. The space observatory Gaia performs a deep optical survey with exquisite astrometric accuracy, for both position and proper motions, that can be used to study natal kicks in NS HMXBs. We aim to survey the observed Galactic NS HMXB population and to quantify the magnitude of the kick imparted onto their NSs, and to highlight any possible differences arising in between the various HMXB types. We perform a census of Galactic NS HMXBs and cross-match existing detections in X-rays, optical and infrared with the Gaia Early Data Release 3 database. We retrieve their parallaxes, proper motions, and radial velocities (when available), and perform a selection based on the quality of the parallax measurement. We then compute their peculiar velocities with respect to the rotating reference frame of the Milky Way, and including their respective masses and periods, we estimate their kick velocities through Markov Chain Monte Carlo simulations of the orbit undergoing a supernova event. We infer the posterior kick distributions of 35 NS HMXBs. After an inconclusive attempt at characterising the kick distributions with Maxwellian statistics, we find that the observed NS kicks are best reproduced by a Gamma distribution of mean 116^+18^_-15_km/s. We note that supergiant systems tend to have higher kick velocities than Be High-Mass X-ray Binaries. The peculiar velocity versus non-degenerate companion mass plane hints at a similar trend, supergiant systems having a higher peculiar velocity independently of their companion mass.