We modeled gravity data to explore Mercury’s internal structure and show the presence of crustal heterogeneities in density. We first evaluated the lithospheric flexure occurring in the spherical harmonic degree range 5–80, according to the flexural isostatic response curve. We thus estimated a mean elastic lithosphere thickness of about 30 pm 10 km and modeled the crust-mantle interface, which varies from 19 to 42 km depth, according to a flexural compensation model. The isostatic gravity anomalies were then obtained as the residual field with respect to the contributions from topography and lithospheric flexure. Isostatic anomalies are mainly related to density variations in the crust: gravity highs mostly correspond to large-impact basins suggesting intra-crustal magmatic intrusions as the main origin of these anomalies. Isostatic gravity lows prevail, instead, above intercrater plains and may represent the signature of a heavily fractured crust.