Density Functional Theory (DFT) calculations within the Generalized Gradient Approximation (GGA) and the GGA + U approach are carried out to investigate the adsorption of O₂ on anatase (101) surfaces having subsurface oxygen vacancies. Our results show that O₂ adsorption is strongly enhanced at sites close to the subsurface defect, whereas dissociation is unfavorable at all sites. The adsorption is accompanied by the transfer of the defect electrons to O₂-derived electronic states in the anatase surface band gap. Peroxide species (O₂²⁻, O–O = 1.48 Å) are stable when the number of adsorbed O₂ molecules is less or equal the number of defects, whereas superoxide species (O₂⁻, O–O = 1.33 Å) become more favorable at coverages exceeding approximately 1.5 O₂ molecules per oxygen vacancy.