This form of Maxwell's Equations should cover metals dielectrics and most other media. Focus on the following
$$
\begin{align}
\boldsymbol{\nabla}.\mathbf{D}&=\rho_f \\
\boldsymbol{\nabla}\times\mathbf{H}&=\mathbf{J}_f+\partial_t\mathbf{D}
\end{align}
$$
Now take divergence of the second equation, and bear in mind that divergence of a curl is zero:
$$
\boldsymbol{\nabla}.\boldsymbol{\nabla}\times\mathbf{H}=\boldsymbol{\nabla}.\mathbf{J}_f+\boldsymbol{\nabla}.\partial_t \mathbf{D}=0
$$
You can then swap and substitute $\boldsymbol{\nabla}.\partial_t \mathbf{D}=\partial_t \boldsymbol{\nabla}.\mathbf{D}=\partial_t \rho_f$
So the continuity equation is a direct consequence of Maxwell's equations. Verifying it is akin to verifying Maxwell's equations. No need to mess around with thing wires and loops, the above prescription applies to most media and even mixed domains, i.e. with different media.
The above analysis tells you that you cannot talk about charge density and current density separately.
I am not sure what proof you are expecting for infinitely thin wire. There I would simply state current density as:
$$
\mathbf{J}\left(\mathbf{r}\right)=\alpha\int ds'\: \boldsymbol{\mathcal{\dot{R}}}\left(s'\right)\delta^{\left(3\right)}\left(\mathbf{r}-\boldsymbol{\mathcal{R}}\left(s'\right)\right)
$$
Where $\boldsymbol{\mathcal{R}}\left(s\right)$ is the trajectory of your wire, parametrized by arc-length $s$. Assuming $\boldsymbol{\mathcal{\dot{R}}}.\boldsymbol{\mathcal{\dot{R}}}=\mathcal{\dot{R}}^2=const$
Value of $\alpha$ follows from:
$$
\int_S{d^2}r\: \mathbf{\hat{n}}.\mathbf{J}=I\Delta\left(s\in S\right)\cos\theta_\mathcal{\dot{R}}=\alpha\mathcal{\dot{R}}^2\,\Delta\left(s\in S\right)
$$
Where the integral above is over the surface area $S$ of the dot-product of the current density with the surface normal $\mathcal{\hat{n}}$. Quantity $\Delta\left(s\in S\right)=1$ is the wire only goes through the surface area $S$ once. $\cos\theta_{\mathcal{\dot{R}}}$ is the cosine of the angle between the wire and the surface normal. $I$ is the current in the wire. Equation above only makes sense if wire goes through $S$ once or not at all.
From this you can then extract the charge density by taking divergence. I think, you will find that once you take divergence the derivative (of the divergence) and the integral ($\int ds'$) will cancel out, and you will get zero charge density for closed-loop wires
ADDENDUM
Lets see how certain integrals/derivatives will transform
Let us restrict our attention to a region, $\left(x,y,z\right)\in\Omega$ and $s\in\left(s_0,s_1\right)$, where $z$ has a one-to-one relationship with arc-length $s$. More specifically, $\mathcal{R}_z\left(s\right)$ gives the z-coordinate of the arc at all points. Within the region we are considering, let $\mathcal{R}_z^{-1}\left(z\right)$ be defined.
Irrespective of where $z$ is the component of the arc-tangent (of the wire) in the direction of $z$ is:
$$
\frac{d\mathcal{R}_z}{ds}=\mathbf{\hat{z}}.\frac{d\boldsymbol{\mathcal{R}}}{ds}=\mathbf{\hat{z}}.\boldsymbol{\dot{\mathcal{R}}}
$$
Next note that within $\Omega$ for any function $f=f\left(x,y,z\right)$:
$$
\begin{align}
\int^b_a dz\, f\left(x,y,z\right)&=\int^{\mathcal{R}_z^{-1}\left(b\right)}_{\mathcal{R}_z^{-1}\left(a\right)} \frac{dz}{ds} ds f\left(x,y,\mathcal{R}_z\left(s\right)\right)=\int^{\mathcal{R}_z^{-1}\left(b\right)}_{\mathcal{R}_z^{-1}\left(a\right)} \frac{d\mathcal{R}_z}{ds} ds f\left(x,y,\mathcal{R}_z\left(s\right)\right)\\
&=\mathbf{\hat{z}}.\int^{\mathcal{R}_z^{-1}\left(b\right)}_{\mathcal{R}_z^{-1}\left(a\right)} \boldsymbol{\mathcal{\dot{R}}} f\left(x,y,\mathcal{R}_z\left(s\right)\right)ds
\end{align}
$$
Where I simply replaced the integration variable $z\to s$ where $z=\mathcal{R}_z\left(s\right)$. The same trick will work in volume integrals within $\Omega$. The the transform would be:
$$
\begin{align}
x &\to x \\
y&\to y \\
z &\to s \\
\end{align}
$$
With the Jacobian $\left|\frac{\partial\left(x,y,z\right)}{\partial\left(x,y,s\right)}\right|=\mathbf{\hat{z}}.\boldsymbol{\mathcal{\dot{R}}}$
It then follows that (within $\Omega$):
$$
\begin{align}
\mathbf{\hat{z}}.\mathbf{J}\left(\mathbf{r}\right)&=\alpha\mathbf{\hat{z}}.\int ds'\: \boldsymbol{\mathcal{\dot{R}}}\left(s'\right)\delta^{\left(3\right)}\left(\mathbf{r}-\boldsymbol{\mathcal{R}}\left(s'\right)\right)=\alpha\int dz'\: \delta^{\left(3\right)}\left(\left(\begin{array}\\x\\y\\z\end{array}\right)-\left(\begin{array}\\\mathcal{R}_x\left(\mathcal{R}_z^{-1}\left(z'\right)\right)\\\mathcal{R}_y\left(\mathcal{R}_z^{-1}\left(z'\right)\right)\\z'\end{array}\right)\right)=\\
&=\alpha \:\delta\left(x-\mathcal{R}_x\left(\mathcal{R}_z^{-1}\left(z\right)\right)\right)\:\delta\left(y-\mathcal{R}_y\left(\mathcal{R}_z^{-1}\left(z\right)\right)\right)
\end{align}
$$
From here it should be relatively easy to derive:
$$
\int_{x_0}^{x_1} dx \int_{y_0}^{y_1} dy \int_{-l/2}^{l/2} dz\: \mathbf{\hat{z}}.\mathbf{J}\left(\mathbf{r}\right)=l\,\alpha
$$
Assuming that $\left(x_0,\,x_1\right)\times\left(y_0,y_1\right)\times\left(-l/2,l/2\right)$ contain the wire and are within $\Omega$. I think this is what you were after.
Another interesting thing to try is:
$$
\begin{align}
\boldsymbol{\nabla}.\mathbf{J}&=\int^{s_b}_{s_a} ds'\: \boldsymbol{\mathcal{\dot{R}}}\left(s'\right).\boldsymbol{\nabla}\delta^{\left(3\right)}\left(\mathbf{r}-\boldsymbol{\mathcal{R}}\left(s'\right)\right)=\int^{s_b}_{s_a} ds'\: \frac{d\boldsymbol{\mathcal{R}}}{ds'}.\boldsymbol{\nabla}\delta^{\left(3\right)}\left(\mathbf{r}-\boldsymbol{\mathcal{R}}\left(s'\right)\right)
\\
&=-\int^{s_b}_{s_a} ds'\: \frac{d}{ds'}.\delta^{\left(3\right)}\left(\mathbf{r}-\boldsymbol{\mathcal{R}}\left(s'\right)\right)=\delta^{\left(3\right)}\left(\mathbf{r}-\boldsymbol{\mathcal{R}}\left(s_a\right)\right)-\delta^{\left(3\right)}\left(\mathbf{r}-\boldsymbol{\mathcal{R}}\left(s_b\right)\right)
\end{align}
$$
Clearly, if $s_a=s_b$ as would be in the case of a closed loop, the divergence of the current density would vanish.
Best Answer
We have, for a point charge $q$ at position $\vec r(t)$: $$\rho(\vec x, t) = q\delta^3(\vec x - \vec r(t))$$ $$\vec J(\vec x, t) = q \frac{d\vec r}{dt}\delta^3(\vec x - \vec r(t))$$ Let us for now work without worrying about what the derivative (more precisely, gradient) of the delta function actually is. We will also enforce the convention that $\vec\nabla$ denotes differentiation with respect to $\vec x$ unless specified otherwise.
We will first evaluate the rate of change of the charge density: $$\frac{\partial \rho}{\partial t} = \frac{\partial}{\partial t}\left[q\delta^3(\vec x - \vec r(t))\right]$$ $$ = q\frac{d\vec{r}}{dt}\cdot\vec\nabla_r\delta^3(\vec x - \vec r(t))$$ as all the time dependence is contained (through $\vec r(t)$) in the delta function. Now, for $\delta(y-z)$, we have the result that $$\frac{\partial}{\partial y} \delta(y-z) = -\frac{\partial}{\partial z} \delta(y-z)$$ We can use this here to get: $$\frac{\partial \rho}{\partial t} = -q\frac{d\vec{r}}{dt}\cdot\vec\nabla\delta^3(\vec x - \vec r(t))$$
The divergence of the current density is: $$\vec\nabla\cdot\vec J = \vec\nabla\cdot\left[q\frac{d\vec{r}}{dt}\delta^3(\vec x - \vec r(t))\right]$$ $$ = q\frac{d\vec{r}}{dt}\cdot\vec\nabla\delta^3(\vec x - \vec r(t))$$ where the field-like dependence on position is once again contained only in the delta function.
Now, whenever $\vec x \ne \vec r(t)$, we know that any derivative of the delta function (which is uniformly zero in this region) must be zero and the continuity equation trivially holds for this case. However, at $\vec x = \vec r(t)$, it still holds because of the coefficients of the derivative of the delta function in the above expressions (Note: we can treat the delta function as the limiting case of a well-behaved, but sharply peaking function, and it is this consideration that underlies this argument).
Therefore, the equation of continuity holds for all $\vec x$ (and $t$): $$\frac{\partial \rho}{\partial t} + \vec\nabla\cdot\vec J= 0$$