Generally speaking is not possible to assert if the energy is conserved in General Relativity (GR). There are several subtle points about the definition of the energy of the gravitational field and how this could introduce a concept of total energy (including gravitational energy), however, here I will discuss only the energy of the matter content.
In some cases one can prove that the total energy of the matter content is indeed conserved. The total energy momentum tensor (EMT) $T_{\mu\nu}$ must satisfy $$\nabla_\mu{}T^\mu{}_\nu = 0.$$ These conditions come from the Bianchi identities along with the Einstein's equations. Given some timelike vector field $v^\mu$ it is possible to define the energy momentum flow through the foliation defined by $v^\mu$ as $P^\mu = T^{\mu\nu}v_\nu$ (not all timelike vector fields define a global foliation, however, I will ignore this point here).
The vector field $P^\mu$ has its divergence given by $$\nabla_\mu{}P^\mu = \nabla_\mu{}T^{\mu\nu}v_\nu + T^{\mu\nu}\nabla_\mu{}v_\nu = T^{\mu\nu}\nabla_\mu{}v_\nu = T^{\mu\nu}\nabla_{(\mu}{}v_{\nu)},$$ where the parenthesis represent the symmetric part and the last equality comes from the symmetric property of the EMT. If for some reason $\nabla_\mu{}P^\mu = 0$, then the Stoke's theorem (plus some conditions on the manifold or on the EMT at infinity) guarantees that the total energy is conserved, i.e., $$\left.\int\mathrm{d}^3x\sqrt{\gamma}P^\nu{}v_\nu\right\vert_{t_1} = \left.\int\mathrm{d}^3x\sqrt{\gamma}P^\nu{}v_\nu\right\vert_{t_2},$$ where $\gamma$ is the determinant of the metric projected on the spatial hypersurface defined by $v^\nu$ and $t_1, t_2$ are two labels defining two different hypersurfaces.
If $v_\mu$ is a Killing vector field it satisfy $$\mathcal{L}_v g_{\mu\nu} = n^\alpha\nabla_\alpha{}g_{\mu\nu} + 2\nabla_{(\mu}v_{\nu)} = 2\nabla_{(\mu}v_{\nu)} = 0,$$ where we used the covariant derivative compatible with $g_{\mu\nu}$. This shows that, if there is a timelike Killing field then the total energy is conserved, however, the converse is not true, i.e., the following statement is not true: if the total energy momentum is conserved then there is a timelike Killing field.
In our current cosmological model, the universe (at zero order) is described by a Friedmann-Lamaître-Robertson-Walker (FLRW) metric which do not posses a timelike Killing vector. This is why the energy of a radiation-like fluid is not conserved, from the direct calculation of the EMT divergence in a FLRW model we have $$\dot{\rho} + 3H(\rho+p) = 0,$$ where $H=\dot{a}/a$ is the Hubble function, $a$ the scale factor, $\dot{f} = v^\mu\partial_\mu{}f$ the time derivative of a scalar function $f$, the EMT is given by $T_{\mu\nu} = \rho{}v_\mu{}v_\nu + p\gamma_{\mu\nu}$, $v_\mu$ is the field represent the fluid flow, $\gamma_{\mu\nu} = g_{\mu\nu} + v_{\mu}v_{\nu}$ is the spatial projector, $\rho$ the energy density in this frame and $p$ the isotropic pressure also in this frame. For a constant equation of state ($w = p/\rho$) we have $$\rho = \rho_0\left(\frac{a_0}{a}\right)^{3(1+w)},$$ where $a_0$ and $\rho_0$ are the scale factor and energy density calculated in a spatial section defined by $t_0$.
From the direct calculation of the total energy (in this frame) we have $$\int\mathrm{d}^3x\sqrt{\gamma}P^\nu{}v_\nu = \int\mathrm{d}^3x\sqrt{\gamma}\rho = \int\mathrm{d}^3x\sqrt{\gamma_0}\rho_0\left(\frac{a_0}{a}\right)^{3w},$$ where we used that $\dot{\sqrt{\gamma}} = 3H\sqrt{\gamma}$ and, therefore, $\sqrt{\gamma} = \sqrt{\gamma_0}(a_0/a)^3$ (this is true in FLRW, in general $3H$ is substituted by the expansion factor $\Theta \equiv \nabla_\mu{}v^\mu$). This shows that, for radiation ($w=1/3$) the energy decreases with $\propto a^{-1}$ when the universe expands. Note also that the presence of dark energy $w<-1/3$ make the total energy increases, e.g., the cosmological constant has $w=-1$ then the energy goes like $a^3$.
For the special case of dust $w=0$ the total energy is conserved. This is an example of what I said before, we can have total energy conservation without a Killing field, in this case this happens because $T_{\mu\nu} = \rho{}v_\mu{}v_\nu$ is orthogonal to $\nabla_\mu{}v_\nu = \mathcal{K}_{\mu\nu} = H\gamma_{\mu\nu}$, where $\mathcal{K}_{\mu\nu}$ is the extrinsic curvature that in FLRW is proportional to $\gamma_{\mu\nu}$.
Finally, we can only have thermodynamic equilibrium when we have a timelike Killing field with the exception that for radiation we just need a conformal Killing field to achieve equilibrium (see "Kinetic theory in the expanding universe" Bernstein 1988). In a FLRW universe we have a timelike conformal Killing field and thats why we have a well defined temperature for radiation, using the Bose-Einstein distribution (assuming kinetic equilibrium) we obtain that $T \propto a^{-1}$, thats why, in the thermodynamical viewpoint, the total energy is not conserved, the temperature drops when the universe expands.
