Mid-lithosphere discontinuities are seismic interfaces likely located within the lithospheric mantle of stable cratons, which typically represent velocities decreasing with depth. The origins of these interfaces are poorly understood due to the difficulties in both characterizing them seismically and reconciling the observations with thermal-chemical models of cratons. Metasomatism of the cratonic lithosphere has been reported by numerous geochemical and petrological studies worldwide, yet its seismic signature remains elusive. Here, we identify two distinct mid-lithosphere discontinuities at ~89 and ~115 km depth beneath the eastern Wyoming craton and the southwestern Superior craton by analyzing seismic data recorded by two longstanding stations. Our waveform modeling shows that the shallow and deep interfaces represent isotropic velocity drops of 2–9% and 3–10%, respectively, depending on the contributions from changes in radial anisotropy and density. By building a thermal-chemical model including the regional xenolith thermobarometry constraints and the experimental phase-equilibrium data of mantle metasomatism, we show that the shallow interface probably represents the metasomatic front, below which hydrous minerals such as amphibole and phlogopite are present, whereas the deep interface may be caused by the onset of carbonated partial melting. The hydrous minerals and melts are products of mantle metasomatism, with CO2-H2O-rich siliceous melt as a probable metasomatic reagent. Our results suggest that mantle metasomatism is probably an important cause of mid-lithosphere discontinuities worldwide, especially near craton boundaries, where the mantle lithosphere may be intensely metasomatized by fluids and melts released by subducting slabs.