The presently available staggered-grid finite-difference (SGFD) schemes for the 3D first-order elastic-wave equation can only achieve high-order spatial accuracy, but they exhibit second-order temporal accuracy. Therefore, the commonly used SGFD methods may suffer from visible temporal dispersion and even instability when relatively large time steps are involved. To increase the temporal accuracy and stability, we have developed a novel time-space-domain high-order SGFD stencil, characterized by ()th-order spatial and ()th-order temporal accuracies (), to numerically solve the 3D first-order elastic-wave equation. The core idea of this new stencil is to use a double-pyramid stencil with an operator length parameter together with the conventional second-order SGFD to approximate the temporal derivatives. At the same time, the spatial derivatives are discretized by the orthogonality stencil with an operator length parameter . We derive the time-space-domain dispersion relation of this new stencil and determine finite-difference (FD) coefficients using the Taylor-series expansion. In addition, we further optimize the spatial FD coefficients by using a least-squares (LS) algorithm to minimize the time-space-domain dispersion relation. To create accurate and reasonable P-, S-, and converted wavefields, we introduce the 3D wavefield-separation technique into our temporal high-order SGFD schemes. The decoupled P- and S-wavefields are extrapolated by using the P- and S-wave dispersion-relation-based FD coefficients, respectively. Moreover, we design an adaptive variable-length operator scheme, including operators and , to reduce the extra computational cost arising from adopting this new stencil. Dispersion and stability analyses indicate that our new methods have higher accuracy and better stability than the conventional ones. Using several 3D modeling examples, we demonstrate that our SGFD schemes can yield greater temporal accuracy on the premise of guaranteeing high-order spatial accuracy. Through effectively combining our new stencil, LS-based optimization, large time step, variable-length operator, and graphic processing unit, the computational efficiency can be significantly improved for the 3D case.