Best representation of a proton in non-static spacetimes The classic way to represent a proton in flat spacetime is to use a simple point charge,
$$J = (e \delta(r), 0)$$
But if we're using an arbitrary metric (let's assume still globally hyperbolic orientable to avoid any issues), then the local conservation of charge is expressed via
\begin{eqnarray}
\nabla_\mu J^\mu &=& \partial_\mu J^\mu + \frac{\partial_\mu \sqrt{g}}{\sqrt{g}} J^\mu\\
&=& \dot{\rho} + \vec{\nabla} \cdot \vec{J} + \frac{\partial_t \sqrt{g}}{\sqrt{g}} \rho + \frac{\vec{\nabla} \sqrt{g}}{\sqrt{g}} \cdot \vec{J}\\
&=& 0
\end{eqnarray}
Obviously if our spacetime isn't static, we get that 
\begin{eqnarray}
\frac{d}{dt} Q = - \int_\Sigma \left[ \vec{\nabla} \cdot \vec{J} + \frac{\partial_t \sqrt{g}}{\sqrt{g}} \rho + \frac{\vec{\nabla} \sqrt{g}}{\sqrt{g}} \cdot \vec{J} \right] \sqrt{q} d^3x
\end{eqnarray}
With $q$ the metric on the Cauchy surface $\Sigma$. Given enough conditions on the metric, we could probably get back conservation of total charge, using Stokes theorem, but for a non-static one, it seems a forlorn hope.
What would then be the most appropriate for the notion of a point charge then? Do we simply pick $\rho = e \delta(r)$, and have the total charge vary, or do we pick a charge density such that the total charge is conserved? I'm not sure which one would approximate the behaviour of a real charge best here. This may also make it depend on the choice of foliation so I'm not sure if the notion of total charge is even that decent in this case.
 A: The metric isn't a separate thing from the existence of this charge. In general, you need a consistent solution to the field equations that has a charge in it. You aren't going to get that if you try to describe the proton as a point charge. For the reasons why this doesn't work, see this question: what is the difference between a blackhole and a point particle
A spacetime describing a proton would contain a charged matter field localized to a region of radius ~1 fm. The stress-energy tensor would contain a contribution from the electric field of the proton, and also a contribution from the force holding the proton together. IIRC there are various more or less messy ways of approximating this type of solution in closed form. You take some equation of state for a sphere, match boundary conditions, etc. I've seen it done for a self-gravitating body like the earth.
When you put such a thing on a background, you're going to run into the problem that GR is nonlinear, so you can't in general just put together two solutions. It's not like E&M, where you can form a solution to Maxwell's equations for a point charge, and another solution representing some external field, and just add them to find a new solution. But this is not a big problem in the approximation that the scale of the background curvature is large compared to the size of the proton.

This may also make it depend on the choice of foliation so I'm not sure if the notion of total charge is even that decent in this case.

I don't think this is a big issue. GR can't formulate global conservation laws for vector quantities like energy-momentum, because parallel transport of vectors is path-dependent, and therefore you can't just say, "let's add up all the energy-momentum on this Cauchy surface." But charge is a scalar. There is a good discussion of this kind of thing in Misner, Thorne, and Wheeler, p. 457. They discuss the example of trying to define the total charge of the universe when the topology is closed -- this example ends up not being well defined.
I don't think there is any huge problem in general with modeling E&M in curved spacetime. You just write down Maxwell's equations in covariant form. But I think there are difficult technical issues in some cases. E.g., people are currently working on trying to find violations of cosmic censorship when you dribble charged particles onto a black hole, and there are difficulties in getting a reliable calculation of the back-reaction on the particles.
A: For a point charge $e$  moving along the worldline $x^\mu_e (s)$ the current density is
$$
J^\mu (x) = e \int u^\mu \frac{\delta^{4}(x-x_e(s))}{\sqrt{-g}} ds,
$$
where the 4D delta-function is defined via $\int \delta^{4}(x)d^4x=1$, i.e. non-covariantly.
For a coordinate system where a point charge has fixed spatial coordinates $\mathbf{x}_e$ and timelike coordinate $x^0$ coincides with its proper time we would have
$$
J^{\mu}=\left(e \frac{\delta^3 (\mathbf{x}-\mathbf{x}_e)}{\sqrt{-g}},0\right).
$$
