Taking a GR perspective, the goal would be to solve the coupled Einstein-Maxwell equations for cosmological boundary conditions.
In the usual Einstein-only case (no net charge density), the Einstein equations, assuming homogeneity and isotropy, reduce to the Friedmann equations. At a high level, the way this derivation works is to first note that the symmetries of the problem mean that (ignoring spatial curvature since it isn't really relevant for this question) we can choose coordinates where the metric has the form
$$
ds^2 = -dt^2 + a^2(t) \delta_{ij} d x^i d x^j
$$
where $\delta_{ij}$ is the Kronecker delta, and the $i,j$ indices run over $1,2,3$. Furthermore, we can use the symmetries to assume the stress energy tensor has the form
$$
T^{0}_{\ \ 0} = -\rho, \ \ T^i_{\ \ j} = p \delta^i_j, \ \ T^0_{\ \ i}= 0
$$
I want to make the the point that both of the above decompositions rely on the symmetry in two ways. First, we can do a split of space and time such that there is no mixing between them -- formally, we can choose a gauge where the shift vector is zero, less formally, there are no $0i$ components for the tensors. Second, on a spatial slice, the spatial components of the tensors are proportional to $\delta^i_j$, because that is the only tensorial object that is invariant under rotations on the spatial slice. Given the above ansatz, we then plug everything into Einstein's equations, and find they reduce to the Friedmann equations, which relate $a(t)$ with $\rho$ and $p$.
Now if we have a charge distribution, we will also need to include the current $J^\mu$ and electromagnetic stress tensor $F_{\mu\nu}$. There is no symmetry reason we can't have $J^0\neq 0$ and $J^i=0$. However, we cannot do the same split for $F_{\mu\nu}$. Because $F_{\mu\nu}$ is antisymmetric, $F_{00}=0$, and if we require $F_{ij}= f(t) \delta_{ij}$ for some $f(t)$, it must be that $f(t)=0$ since $F_{ij}$ is antisymmetric but $\delta_{ij}$ is symmetric. On top of that, there are no two index anti-symmetric tensors that are invariant under rotations that we can use to be proportional to $F_{\mu\nu}$.
As a result, starting from the cosmological assumptions of homogeneity and isotropy, we are stuck if we try to find a non-trivial background solution involving $F_{\mu\nu}$. There isn't an appropriate tensor structure that can support a non-vanishing $F_{\mu\nu}$. Less formally, the symmetries of the problem are highly constraining and force $F_{\mu\nu}=0$.
You could get around this by breaking the strong requirements of exact homogeneity and isotropy. You could imagine that homogeneity is broken, and there are patches of positive and negative charges, for example. So I wouldn't say the result is that the Universe can't have a net charge. Just that there's no way to have a non-trivial electromagnetic field in a Universe that is exactly homogenous and isotropic.