Reservoir depletion can induce substantial changes in the stress state of the rock. The coupled interaction between the pore fluid pressure and rock stress will then alter the reservoir permeability, which in turn reversely affects the productivity index of the production well. A new nonlinear analytical solution is developed for the drawdown-dependent productivity index of reservoirs under steady-state flow. Biot's theory of poroelasticity is used to derive the depletion-induced changes in the reservoir rock porosity and permeability. The well-known Mindlin's solution for a Nucleus of Strain in a semi-infinite elastic medium is applied as Green's function and integrated over the depleted volume of reservoir rock to obtain the 3D distribution of stress and volumetric strain distributions. The fluid transport equation is nonlinearly coupled to the solid mechanics solution via the stress-dependent permeability coefficients. A perturbation technique is applied to mathematically treat the described nonlinearity to solve for the coupled equations of pore fluid flow and rock stress under steady-state flow. The good match between the obtained analytical approximations for productivity index and the numerical solutions verifies the correctness and robustness of the proposed model. Results indicate and confirm the expected strong dependency of the well productivity index to the drawdown magnitude as well as the poroelastic constitutive parameters of the reservoir rock, with the highest sensitivity to drained bulk modulus, followed by the reservoir depth and solid-grain modulus. The lowest PI sensitivity is to the pore fluid modulus and Poisson's ratio. The resulting productivity index is found out to be drawdown-dependent, which can render values substantially different than the productivity index estimate from the conventional flow-only analysis. The presented estimates for the related nonlinear productivity index can be readily used by the practicing engineers.