Direct numerical simulations (DNS) were carried out to investigate turbulent heat transfer in a channel flow of homogenous polymer solutions described by the Finitely Extensible Nonlinear Elastic-Peterlin (FENE-P) constitutive model at intermediate and high Prandtl numbers (Pr = 1.25 and 5). Time averaged statistics of temperature fluctuations, turbulent heat fluxes, thermal turbulent diffusivity, and of budget terms of the temperature variance are reported and compared with those of the Newtonian fluid cases at the same Prandtl and Reynolds numbers. Moreover, twenty one sets of DNS data of fluid flow are utilized to improve existing k-epsilon-(nu(2)) over bar -f models for FENE-P fluids to deal with turbulent flow of dilute polymer solutions up to the high drag reduction regime; specifically the dependency of closures on the wall friction velocity is removed. Furthermore, five sets of recent DNS data of fluid flow and heat transfer of FENE-P fluids were used to develop the first RANS model capable of predicting the heat transfer rates in viscoelastic turbulent flows. In this model, an existing closure for calculating the turbulent Prandtl number for Newtonian fluids is extended to deal with heat transfer in turbulent viscoelastic fluids. Predicted polymer stresses, velocity profiles, mean temperature profiles, and turbulent flow characteristics are all in good agreement with the DNS data, and show improvement over previous RANS models. (C) 2016 Elsevier Ltd. All rights reserved.