Since the discovery of an excess in gamma rays in the direction of M31, its cause has been unclear. Published interpretations focus on dark matter or stellar related origins. Studies of a similar excess in the Milky Way center motivate a correlation of the spatial morphology of the signal with the distribution of stellar mass in M31. However, a robust determination of the best theory for the observed excess emission is challenging due to uncertainties in the astrophysical gamma-ray foreground model. We perform a spectro-morphological analysis of the M31 gamma-ray excess using state-of-the-art templates for the distribution of stellar mass in M31 and novel astrophysical foreground models for its sky region. We construct maps for the old stellar populations of M31 based on data from the PAndAS survey and carefully remove the foreground stars. We also produce improved astrophysical foreground models via novel image inpainting techniques based on machine learning methods. Our stellar maps, mimicking the location of a population of millisecond pulsars in the bulge of M31, reach a 5.4σ significance, making them as strongly favoured as the simple phenomenological models usually considered in the literature, e.g. disk-like templates. This detection is robust to generous variations of the astrophysical foreground model. Once the stellar templates are included in the astrophysical model, we show that the dark matter annihilation interpretation of the signal is unwarranted. We demonstrate that about one million unresolved millisecond pulsars naturally explain the observed gamma-ray luminosity per stellar mass, energy spectrum, and stellar bulge-to-disk flux ratio.