Turbulent heat transfer to CO2 at supercritical pressure flowing in heated vertical tubes is investigated using direct numerical simulation at the inlet Reynolds number Re0=5400, which is based on inlet bulk velocity and tube diameter. Temperature range within the flow field covers the pseudocritical region, where very significant fluid property variations are involved. Both upward and downward flows are considered. The wall temperature distribution shows well-known heat transfer deterioration characterized by the localized peak in upward flows, while no such anomaly is observed in downward flows. The deterioration occurs at the region where turbulence is attenuated significantly, and is followed by the enhancement with restoration of turbulence caused by complicated interactions with a buoyancy effect. Further investigation of turbulence statistics indicates that ρux″ux″¯, ρux″ur″¯, and ρux″h″¯ are significantly affected by their respective buoyancy production terms due to ρ′ux′¯, ρ′ur′¯, and ρ′h′¯ which are proven to be significant in vertical supercritical flows. Combined with the deformation of mean velocity profile into an M-shaped one in upward flow, ρ′ux′¯ becomes negatively correlated from a certain downstream region so that ρux″h″¯ undergoes a very complicated transition changing both in sign and magnitude, causing severe impairment of heat transfer in upward supercritical flows.