Depth-dependent density perturbations are determined inside the elastic part of the Martian lithosphere (assumed to be 300 km thick) using a thick elastic shell model, such that the deformed model under the loads of surface topography and the density perturbations gives rise to the observed topography and gravity anomalies of Mars specifies by spherical harmonics of degree 3--50. Large positive density perturbations, as high as ~150 kg/m3, exist in the crust that are associated with Hellas, Isidis, Argyre, and Utopia basins, reflecting the mantle uplift in response to the impact excavation. Also, large positive density perturbations exist in the crust beneath Olympus (~350 kg/m3) and Tharsis monts (~200 kg/m3), indicating that these shield volcanoes are not isostatically compensated. On the other hand, negative density perturbations are associated with Alba Patera, Elysium rise, and the highland surrounding Hellas basin. The density perturbations gradually diminish with depth, becoming less than 20 kg/m3 by 250--300 km depth. The resulting stress differences within the crust and upper mantle are dominated by Olympus (up to 100 MPa), whereas other shield volcanoes create lesser stress differences. We estimate the brittle-ductile transition temperature in the lithosphere, adopting the creep law of dry olivine for the rheology of the Martian lithosphere. The brittle-ductile transition temperature beneath Olympus is ~800C throughout the lithosphere, indicating that the upper parts of Mars must have been cold to support this volcano for 1--2 Gyr.