  We provide new insights into the a priori theory for a time-stepping scheme
based on least-squares finite element methods for parabolic first-order
systems. The elliptic part of the problem is of general
reaction-convection-diffusion type. The new ingredient in the analysis is an
elliptic projection operator defined via a non-symmetric bilinear form,
although the main bilinear form corresponding to the least-squares functional
is symmetric. This new operator allows to prove optimal error estimates in the
natural norm associated to the problem and, under additional regularity
assumptions, in the $L^2$ norm. Numerical experiments are presented which
confirm our theoretical findings.
