The crustal structure beneath the Japanese islands, including depth distributions of the Conrad and Moho discontinuities, was estimated using a tomographic inversion of regional body wave arrival times. Depth distributions of the bottom of the surface layer, the Conrad, and the Moho were modeled with two-dimensional B spline functions, while velocity distributions in layers were expressed by three-dimensional B spline functions. The depth of the discontinuities and the velocity in the layers were estimated simultaneously by the least squares method. The velocity structure was sequentially estimated from shallower parts to deeper parts to avoid correlation between them. This sequential analysis provided improved depth resolution. The deepest region of the Moho discontinuity was located in central Honshu, reaching about 40 km. The Moho discontinuity was generally deep in the central part of the islands, whereas it was relatively shallow in the Kanto, southwestern Chubu, and Chugoku districts in Honshu and in northern Kyushu. Some of the shallow Moho regions would be related to graben formation due to tensile tectonic stress since the Miocene. The results were compared with those of seismic refraction surveys and receiver function analyses, and it was found that the obtained model was consistent with many of these studies.