By means of high-resolution cosmological hydrodynamical simulations of Milky Way-like disc galaxies, we conduct an analysis of the associated stellar metallicity distribution functions (MDFs). After undertaking a kinematic decomposition of each simulation into spheroid and disc sub-components, we compare the predicted MDFs to those observed in the solar neighbourhood and the Galactic bulge. The effects of the star formation density threshold are visible in the star formation histories, which show a modulation in their behaviour driven by the threshold. The derived MDFs show median metallicities lower by 0.2-0.3 dex than the MDF observed locally in the disc and in the Galactic bulge. Possible reasons for this apparent discrepancy include the use of low stellar yields and/or centrally-concentrated star formation. The dispersions are larger than the one of the observed MDF; this could be due to simulated discs being kinematically hotter relative to the Milky Way. The fraction of low metallicity stars is largely overestimated, visible from the more negatively skewed MDF with respect to the observational sample. For our fiducial Milky Way analog, we study the metallicity distribution of the stars born in situ relative to those formed via accretion (from disrupted satellites), and demonstrate that this low-metallicity tail to the MDF is populated primarily by accreted stars. Enhanced supernova and stellar radiation energy feedback to the surrounding interstellar media of these pre-disrupted satellites is suggested as an important regulator of the MDF skewness.