Okay, so we are using coordinates $(t,r,\theta,\phi)$, where $t$ and $r$ have dimensions of time and length, and $\theta $ and $\phi$ are dimensionless. We see that $g_{\mu\nu} \mathrm dx^\mu \mathrm dx^\nu$ has consistent units of length squared, as required.
The Stress Energy Momentum Tensor of a perfect fluid from this article is $$T^{\mu\nu}=\left(\rho+\frac{p}{c^2}\right)u^\mu u^\nu-p g^{\mu\nu}$$ with $u=\{1,0,0,0\}$.
This is a mistake. From the article you quote,
[...] where $u$ is the 4-velocity of the fluid and $\eta=\mathrm{diag}(1,-1,-1,-1)$ is the metric tensor of Minkowski spacetime.
This indicates that the expression given is derived in coordinates where $x^0 = ct$ rather than just $t$. As a result, the $00$-component of the metric should be $1$, not $c^2$, and $u=(c,0,0,0)$, not $(1,0,0,0)$. From there, we would have
$$T^{\mu\nu}= \begin{bmatrix}\rho c^2 &0&0&0\\0&\frac{p(1-kr^2)}{a^2}&0&0\\0&0&\frac{p}{a^2 r^2}&0\\0&0&0&\frac{p}{a^2r^2\sin^2(\theta)}\end{bmatrix}$$ $$T_\mu^{\ \ \nu}= \begin{bmatrix}\rho c^2 &0&0&0\\0&-p&0&0\\0&0&-p&0\\0&0&0&-p\end{bmatrix}$$ $$T_{\mu\nu}= \begin{bmatrix}\rho c^2 &0&0&0\\0&\frac{pa^2}{1-kr^2}&0&0\\0&0&pa^2r^2&0\\0&0&0&pa^2r^2\sin^2(\theta)\end{bmatrix}$$
On the other side of the Einstein equations, in coordinates $(ct,r,\theta,\phi)$, the Einstein tensor takes the form
$$G_{00} = 3\left(\frac{\dot a^2+kc^2}{a^2c^2}\right)-\Lambda$$ $$G_{rr} = -\left(\frac{2a \ddot a + \dot a^2 + kc^2}{c^2(1-kr^2)}\right)+\frac{\Lambda a^2}{1-kr^2}$$ $$G_{\theta\theta} = -\frac{2a\ddot a + \dot a^2 + kc^2)}{c^2} r^2 +\Lambda a^2 r^2$$$$G_{\theta\theta} = -\frac{\left(2a\ddot a + \dot a^2 + kc^2\right)}{c^2} r^2 +\Lambda a^2 r^2$$ $$G_{\phi\phi} = -\frac{2a\ddot a + \dot a^2 + kc^2)}{c^2}r^2 \sin^2(\theta)+ \Lambda a^2 r^2\sin^2(\theta)$$$$G_{\phi\phi} = -\frac{\left(2a\ddot a + \dot a^2 + kc^2\right)}{c^2}r^2 \sin^2(\theta)+ \Lambda a^2 r^2\sin^2(\theta)$$
Setting $G_{\mu\nu}= \frac{8\pi G}{c^4} T_{\mu\nu}$, this yields two independent equations: $$3\left(\frac{\dot a^2+kc^2}{a^2c^2}\right)-\Lambda= \frac{8 \pi G}{c^4}(\rho c^2) \iff \frac{\dot a^2}{a^2} + \frac{kc^2}{a^2}= \frac{8\pi G \rho}{3} + \frac{\Lambda c^2}{3} $$ $$-\left(\frac{2a \ddot a + \dot a^2 +k}{c^2(1-kr^2)}\right)+\frac{\Lambda a^2}{1-kr^2}=\frac{8 \pi G}{c^4}\frac{pa^2}{1-kr^2} \iff 2\frac{\ddot a}{a} + \frac{\dot a^2}{a^2} + \frac{kc^2}{a^2} = -\frac{8\pi G}{c^2}p + \Lambda c^2$$
Alternatively, we could use $x^0=t$ as you propose in the beginning, so $g_{00}=c^2$ and $u=(1,0,0,0)$. This amounts to a coordinate change $x=(ct,r,\theta,\phi)\mapsto x'=(t,r,\theta,\phi) \implies \frac{\partial x'^0}{\partial x^0} = 1/c$. Applying the standard tensor transformation rules we would have that $T'^{00}=T^{00}/c^2 = \rho$ and $T'_{00} = T_{00} c^2 = \rho c^4$, which is (if I understand correctly) what you have found.
I've mentioned this ad naseum: natural units are the bane of my existence.
Actually, the problem here is not natural units, but rather choices of coordinate. Whether you use SI or not, you may choose $(ct,r,\theta,\phi)$ or $(t,r,\theta,\phi)$. The most common convention by far is to do the former, so in that sense your choice of the latter will yield differences with most of the established literature.
The key thing to remember is that coordinates - and by extension, the components of tensors - are arbitrary. The value (and even the dimensions) of $T_{00}$ depend on which coordinates you choose, even within the SI system of units. The full expression for the tensor $T$ is $$T = T_{\mu\nu} \mathrm dx^\mu \mathrm dx^\nu$$ and it is this which must have consistent dimensions across all coordinate systems. In SI, it has units of $\mathrm J \mathrm m^{-1} \mathrm s^{-2}$, so the units of $T_{00}$ are $\frac{\mathrm J \mathrm m^{-1} \mathrm s^{-2}}{[x^0]^2}$, where $[a]$ means "units of $a$."
