In general it is hard to assign a total energy to a universe in General Relativity. This is a FAQ <https://math.ucr.edu/home/baez/physics/Relativity/GR/energy_...>. To that we can add that we just don't know enough about our universe -- how big is it? how much matter is in it? -- because we can only see a fraction of the total universe both because a lot of it has expanded away from our view and because our instruments are not yet good enough to pick out features obscured by dust and gas and/or by very low interaction with electromagnetism (doesn't create strong spectral emission or absorption lines).
What we can do is consider the energy-density at a point in a general curved spacetime. That quantity can be subdivided into roughly the part that moves ultra-relativistically (including radiation like photons), that which moves slowly compared to that (including the overwhelming majority of protons and neutrons). From there we can justify a universal average because the universe is (at large scales) homogenous and isotropic: with a wide enough view everything looks the same, and in every direction.
Once we have the average energy-density we can think of how it evolves in an expanding universe. Non-relativistic matter stays mostly clumped in galaxies and those separate, but because we take an average this means that non-relativistic matter's contribution to the energy-density at any point drops with expansion. Ultra-relativistic matter spreads everywhere, so it too dilutes, and additionally it gets colder (seen as e.g. a drop in the energy-density of photons at a point, which corresponds with a drop in a the average photon's frequency or equivalently a stretching of the average photon's wavelength; the latter might lead to an intuitive understanding that an average photon gets smeared over increasing numbers of points, so there is less photon or less photon-energy at any point).
So far this is treating photons like classical particles, or even like discrete quantum particles. In physical cosmology, one would want to consider them as the field content of the electromagnetic field of the Standard Model, which is a relativistic quantum field theory (QFT). Often wanting to do it "properly" runs into tractability problems in calculations, but there are some important features of a QFT compared to a non-relativistic theory where there is some constant number of photons.
The most important part is that last: in a QFT, the count of photons depends on the observer's acceleration. An accelerated observer sees more particles than a freely-falling, zero-angular-momentum, zero-force observer who cannot measure any acceleration. This is the Unruh effect. Both the accelerated observer and the unaccelerated one count correctly, and can calculate what the other observer counts, although each may have different explanations about why its count differs from the other's. The equivalence principle tells us that an observer standing on a gravitating object (like the surface of the Earth) sees most of the effects an observer in deep space accelerated by some propulsive mechanism. The energy-density at a point, in a QFT on curved spacetime, is an observer-dependent quantity.
An accelerating expanding universe has tremendous spacetime curvature, and one result is that the particle-counts obtained by an observer depends on when the particles are counted, all other things being equal. A future observer will count more particles than a past one. (This is also how Hawking radiation works: an observer in the past of an evaporating black hole sees fewer particles than a later observer; those particles are mostly in the region near the black hole; by contrast a late-time cosmological observer will see "extra" or "new" particles at great distances).
So, what can we say about the total photon content? Firstly, it's time-dependent. Secondly, it depends on how the observer is moving against the electromagnetic field in the appropriate QFT. Thirdly, it depends on detailed knowledge of the distribution of photons (practically impossible, so we average). Fourthly it depends on detailed knowledge of the expansion, curvature, and acceleration parameters of the expanding spacetime (in practice that's also hard to measure). Fifthly, it would depend on recovering details of the universe that are most likely inaccessible to us because of the expansion carrying that out of view, but also because the hot dense early days of our universe are electromagnetically opaque, and the hotter denser earlier moments of our universe probably have no photons at all (the start of the "electro-weak epoch").
Consequently, our best modelling is only approximations based on the information we have at any given moment. However, such modelling is amenable to perturbative methods, where we add a small difference "by hand" and trace out the consequences to the modelled universe. Adding a photon here does not require removal of a photon there, and in any rate the model should not keep the total photon count constant for general observers (if you hold it constant for a particular chosen observer or family of observers, it won't be constant for other, different observers). So for "Eulerian" observers freely floating and looking at the expansion flow, adding a small perturbation like an "extra" photon should not make much difference. So instead we might want to treat the new photon as a "Lagrangian" observer, tracing out its peculiar experience of the expanding universe in which it exists.
Taking that latter approach, we find that we can readily introduce a photon such that at its distant past it must have had a truly microscopic wavelength, and thus at the point it occupies it vastly exceeds the average photon-energy density of the universe. When we make a small perturbation we expect only small effects. But in this case, by adding a pretty vanilla photon at a point in the universe's distant future and tracing out the consequence, we find we require a big change -- the energy-density "bump" at the points in our test photon's past gets higher and higher, and it ultimately becomes a significant anisotropy and inhomogeneity that disrupts the justification of our averaging. In extremis, the result would be a part of the universe that stands out as obviously very different from what we see in the sky. Indeed, we can set things up so that we can't even make a prediction from Standard Model physics about the behaviour of the ultra-high-energy early photon that we expect to redshift into our test photon we added in as a small perturbation.
In essence, this is just another example of the difficulties one might encounter when mixing quantum field theory with general relativity. General relativity does not in itself forbid even infinite energy at a point, but the behaviour of the Standard Model (a quantum field theory) is not well understood above about 10^14 GeV/c^2 energies. General relativity requires a constraint imposed upon choices of initial values, which in the case of our test photon it might get from an extension to the Standard Model, or which might be put in "by hand"; the Standard Model needs an ultraviolet completion for it to provide or at least justify the constraint.
In general it is hard to assign a total energy to a universe in General Relativity. This is a FAQ <https://math.ucr.edu/home/baez/physics/Relativity/GR/energy_...>. To that we can add that we just don't know enough about our universe -- how big is it? how much matter is in it? -- because we can only see a fraction of the total universe both because a lot of it has expanded away from our view and because our instruments are not yet good enough to pick out features obscured by dust and gas and/or by very low interaction with electromagnetism (doesn't create strong spectral emission or absorption lines).
What we can do is consider the energy-density at a point in a general curved spacetime. That quantity can be subdivided into roughly the part that moves ultra-relativistically (including radiation like photons), that which moves slowly compared to that (including the overwhelming majority of protons and neutrons). From there we can justify a universal average because the universe is (at large scales) homogenous and isotropic: with a wide enough view everything looks the same, and in every direction.
Once we have the average energy-density we can think of how it evolves in an expanding universe. Non-relativistic matter stays mostly clumped in galaxies and those separate, but because we take an average this means that non-relativistic matter's contribution to the energy-density at any point drops with expansion. Ultra-relativistic matter spreads everywhere, so it too dilutes, and additionally it gets colder (seen as e.g. a drop in the energy-density of photons at a point, which corresponds with a drop in a the average photon's frequency or equivalently a stretching of the average photon's wavelength; the latter might lead to an intuitive understanding that an average photon gets smeared over increasing numbers of points, so there is less photon or less photon-energy at any point).
So far this is treating photons like classical particles, or even like discrete quantum particles. In physical cosmology, one would want to consider them as the field content of the electromagnetic field of the Standard Model, which is a relativistic quantum field theory (QFT). Often wanting to do it "properly" runs into tractability problems in calculations, but there are some important features of a QFT compared to a non-relativistic theory where there is some constant number of photons.
The most important part is that last: in a QFT, the count of photons depends on the observer's acceleration. An accelerated observer sees more particles than a freely-falling, zero-angular-momentum, zero-force observer who cannot measure any acceleration. This is the Unruh effect. Both the accelerated observer and the unaccelerated one count correctly, and can calculate what the other observer counts, although each may have different explanations about why its count differs from the other's. The equivalence principle tells us that an observer standing on a gravitating object (like the surface of the Earth) sees most of the effects an observer in deep space accelerated by some propulsive mechanism. The energy-density at a point, in a QFT on curved spacetime, is an observer-dependent quantity.
An accelerating expanding universe has tremendous spacetime curvature, and one result is that the particle-counts obtained by an observer depends on when the particles are counted, all other things being equal. A future observer will count more particles than a past one. (This is also how Hawking radiation works: an observer in the past of an evaporating black hole sees fewer particles than a later observer; those particles are mostly in the region near the black hole; by contrast a late-time cosmological observer will see "extra" or "new" particles at great distances).
So, what can we say about the total photon content? Firstly, it's time-dependent. Secondly, it depends on how the observer is moving against the electromagnetic field in the appropriate QFT. Thirdly, it depends on detailed knowledge of the distribution of photons (practically impossible, so we average). Fourthly it depends on detailed knowledge of the expansion, curvature, and acceleration parameters of the expanding spacetime (in practice that's also hard to measure). Fifthly, it would depend on recovering details of the universe that are most likely inaccessible to us because of the expansion carrying that out of view, but also because the hot dense early days of our universe are electromagnetically opaque, and the hotter denser earlier moments of our universe probably have no photons at all (the start of the "electro-weak epoch").
Consequently, our best modelling is only approximations based on the information we have at any given moment. However, such modelling is amenable to perturbative methods, where we add a small difference "by hand" and trace out the consequences to the modelled universe. Adding a photon here does not require removal of a photon there, and in any rate the model should not keep the total photon count constant for general observers (if you hold it constant for a particular chosen observer or family of observers, it won't be constant for other, different observers). So for "Eulerian" observers freely floating and looking at the expansion flow, adding a small perturbation like an "extra" photon should not make much difference. So instead we might want to treat the new photon as a "Lagrangian" observer, tracing out its peculiar experience of the expanding universe in which it exists.
Taking that latter approach, we find that we can readily introduce a photon such that at its distant past it must have had a truly microscopic wavelength, and thus at the point it occupies it vastly exceeds the average photon-energy density of the universe. When we make a small perturbation we expect only small effects. But in this case, by adding a pretty vanilla photon at a point in the universe's distant future and tracing out the consequence, we find we require a big change -- the energy-density "bump" at the points in our test photon's past gets higher and higher, and it ultimately becomes a significant anisotropy and inhomogeneity that disrupts the justification of our averaging. In extremis, the result would be a part of the universe that stands out as obviously very different from what we see in the sky. Indeed, we can set things up so that we can't even make a prediction from Standard Model physics about the behaviour of the ultra-high-energy early photon that we expect to redshift into our test photon we added in as a small perturbation.
In essence, this is just another example of the difficulties one might encounter when mixing quantum field theory with general relativity. General relativity does not in itself forbid even infinite energy at a point, but the behaviour of the Standard Model (a quantum field theory) is not well understood above about 10^14 GeV/c^2 energies. General relativity requires a constraint imposed upon choices of initial values, which in the case of our test photon it might get from an extension to the Standard Model, or which might be put in "by hand"; the Standard Model needs an ultraviolet completion for it to provide or at least justify the constraint.