A novel hydro-mechanical coupled analysis is proposed for calculating the pore pressure and stress evolution of fractured vuggy carbonate reservoirs under depletion or injection. Different from the conventional dual- and triple-porosity models in which the vugs are considered as continuous porosity, the vugs are treated as virtual volumes in this study. The fluid exchange at the vug-to-matrix interface is dynamically calculated with time evolution and the pore pressure in the vug is updated through considering both the fluid material balance and the volume change due to the mechanical deformation of the vug. The fluid-mechanical interaction in the rock matrix and natural fractures is calculated based on the framework of Biot's poroelasticity theory. The mechanical and hydraulic interactions between vugs and matrix are maintained and the stress evolution owing to the depletion can be dynamically updated. The results in this study show that the depletion or injection process is greatly affected by the existence of vugs and natural fractures. Natural fractures provide highly conductive channels that lead to the preferential flow paths. The large fluid storage of vugs impedes the pressure diffusion process and brings perturbation to the pressure pattern near the vugs. The response of a vug to the depletion can be divided into three different stages depending on the depletion time. The effects of fluid storage, matrix/fracture permeability, natural fractures and fluid injection are investigated. The proposed hydro-mechanical coupled analysis provides a new modeling method for studying the hydraulic and mechanical behaviors of fractured vuggy carbonate reservoirs.