Extrapolating wavefields and imaging at each depth during three-dimensional recursive wave-equation migration is a time-consuming endeavor. For efficiency, most commercial techniques extrapolate wavefields through thick slabs followed by wavefield interpolation within each thick slab. In this article, we develop this strategy by associating more efficient interpolators with a Fourier-transform-related wavefield extrapolation method. First, we formulate a three-dimensional first-order separation-of-variables screen propagator for large-step wavefield extrapolation, which allows for wide-angle propagations in highly contrasting media. This propagator significantly improves the performance of the split-step Fourier method in dealing with significant lateral heterogeneities at the cost of only one more fast Fourier transform in each thick slab. We then extend the two-dimensional Kirchhoff and Born-Kirchhoff local wavefield interpolators to three-dimensional cases for each slab. The three-dimensional Kirchhoff interpolator is based on the traditional Kirchhoff formula and applies to moderate lateral velocity variations, whereas the three-dimensional Born-Kirchhoff interpolator is derived from the Lippmann-Schwinger integral equation under the Born approximation and is adapted to highly laterally varying media. Numerical examples on the three-dimensional salt model of the Society of Exploration Geophysicists/European Association of Geoscientists demonstrate that three-dimensional first-order separation-of-variables screen propagator Born-Kirchhoff depth migration using thick-slab wavefield extrapolation plus thin-slab interpolation tolerates a considerable depth-step size of up to 72 ms, eventually resulting in an efficiency improvement of nearly 80% without obvious loss of imaging accuracy. Although the proposed three-dimensional interpolators are presented with one-way Fourier extrapolation methods, they can be extended for applications to general migration methods.