Abstract Turbulence has a profound impact on the evolution of gas and dust in protoplanetary disks (PPDs), from driving the collisions and the diffusion of dust grains, to the concentration of pebbles in giant vortices, thus, facilitating planetesimal formation. The vertical shear instability (VSI) is a hydrodynamic mechanism, operating in PPDs if the local rate of thermal relaxation is high enough. Previous studies of the VSI have, however, relied on the assumption of constant cooling rates, or neglected the finite coupling time between the gas particles and the dust grains. Here, we present the results of hydrodynamic simulations of PPDs with the PLUTO code that include a more realistic thermal relaxation prescription, which enables us to study the VSI in the optically thick and optically thin parts of the disk under consideration of the thermal dust-gas coupling. We show the VSI to cause turbulence even in the optically thick inner regions of PPDs in our two- and three-dimensional simulations. The collisional decoupling of dust and gas particles in the upper atmosphere and the correspondingly inefficient thermal relaxation rates lead to the damping of the VSI turbulence. Long-lived anticyclonic vortices form in our three-dimensional simulation. These structures emerge from the turbulence in the VSI-active layer, persist over hundreds of orbits and extend vertically over the whole extent of the turbulent region. We conclude that the VSI leads to turbulence and the formation of long-lived dust traps within ±3 pressure scale heights distance from the disk midplane.