Molecular line observations of starless (prestellar) cores combined with a chemical evolution modeling and radiative transfer calculations are a powerful tool to study the earliest stages of star formation. However, conclusions drawn from such a modeling may noticeably depend on the assumed thermal structure of the cores. The assumption of isothermality, which may work well in chemo-dynamical studies, becomes a critical factor in molecular line formation simulations. We argue that even small temperature variations, which are likely to exist in starless cores, can have a non-negligible effect on the interpretation of molecular line data and derived core properties. In particular, ``chemically pristine isothermal cores (low depletion) can have centrally peaked C$^{18}$O and C$^{34}$S radial intensity profiles, while having ring-like intensity distributions in models with a colder center and/or warmer envelope assuming the same underlying chemical structure. Therefore, derived molecular abundances based on oversimplified thermal models may lead to a mis-interpretation of the line data.