Insights into wave propagation in prestressed porous media are important in geophysical applications, such as monitoring changes in geo-pressure. This can be addressed by poro-acoustoelasticity theory, which extends the classical acoustoelasticity of solids to porous media. The relevant poro-acoustoelasticity equations can be derived from anisotropic poroelasticity equations by replacing the poroelastic stiffness matrix with an acoustoelastic stiffness matrix consisting of second-order and third-order elastic constants. The theory considers the poroelasticity equations to be nonlinear due to the cubic strain-energy function with linear strains under finite-magnitude prestresses. A rotated staggered-grid finite-difference method with an unsplit convolutional perfectly matched layer absorbing boundary is used to solve a first-order velocity-stress formulation of poro-acoustoelasticity equations for elastic wave propagation in prestressed porous media. Numerical solutions are partially verified by computing the velocities of fast P wave, slow P wave, and S wave as a function of hydrostatic prestress and are compared with the exact values. Numerical simulations of wave propagation are carried out for the model of poro-acoustoelastic homogeneous space under three states—prestress confining (hydrostatic), uniaxial, and pure shear—and for the model of two poro-acoustoelastic homogeneous half-spaces in the planar contact under confining (hydrostatic) prestress. The resulting wavefield snapshots show fast P-wave, slow P-wave, and S-wave propagations in poro-acoustoelastic media under loading prestresses, which illustrate that the stress-induced velocity anisotropy is of orthotropy strongly related to the orientation of prestresses. These examples demonstrate the significant impact of prestressing conditions on seismic responses in velocity and anisotropy.