Why does acceleration seem not to be the gradient of gravitational potential? Consider a spherically symmetric distribution of density $\rho(r)$.  We can define the mass enclosed within each radius $r$ using $\frac{dM(r)}{dr} = 4\pi r^2 \rho(r)$, with the condition that $M(r=0) = 0$.
The gravitational acceleration and potential can then be given as,
$$\begin{align}
a_g &= -\frac{GM}{r^2}\\
\Phi_g &= -\frac{GM}{r}.
\end{align}$$
But then the acceleration wouldn't seem to be the gradient of the potential!
$$a_g = - \frac{GM}{r^2} \neq -\nabla \Phi_g = \frac{d}{dr}\left(\frac{GM(r)}{r}\right).$$
What am I missing?

This comes up, for example when writing the equations of stellar structure, or the Lane-Emden equations.  I've seen these written as both,
$\frac{1}{\rho}\frac{dP}{dr} = -\frac{GM}{r^2}$, and as
$\frac{1}{\rho}\frac{dP}{dr} = -\nabla \Phi_g$
for example, the Wikipedia article on Lane-Emden includes both of these expressions.  This seems to suggest that $\Phi_g \neq \frac{GM(r)}{r}$...
 A: Actually, your expression for the potential $\Phi(r)$ is incorrect. The expression
$\Phi(r) = -\frac{GM(r)}{r}$ is only valid outside the sphere.
As an explicit demonstration of its invalidity, note that
$$\underset{r\rightarrow0}{\text{lim}}\,\Phi(r)=\underset{r\rightarrow0}{\text{lim}}\,\left[-\frac{G}{r}\int_0^r4\pi r'^2\rho(r')\,dr'\right]=0$$
assuming that $\rho(r)$ is finite. In particular, this also predicts that $\Phi(0)=0$ for a uniform-density sphere. However, for a uniform density sphere, we actually have $\Phi(0)=-2\pi G\rho R^2$ (using the convention $\Phi(\infty)=0$).
Actual Potential Inside Distribution
The gravitational potential $\Phi_r(a)$ at radius $a$ due to a thin spherical shell of radius $r$ with surface density $\rho(r)$ is
$$\Phi_r(a)=\left\{\begin{aligned}
&-\frac{4\pi G r^2\rho(r)}{a} &&: a>r\\
&-4\pi G \rho(r) r &&: a\leq r
\end{aligned}
\right.$$
and thus the correct expression for the potential due to an arbitrary spherically-symmetric charge distribution $\rho(r)$ becomes
$$\Phi(a)=\int_0^\infty\Phi_r(a)\,dr.$$
If $\rho(r)=0$ for all $r\geq R$, then we can rewrite the integral as
$$\Phi(a)=\int_0^R\Phi_r(a)\,dr=\left\{\begin{aligned}
&\int_0^R-\frac{4\pi G r^2\rho(r)}{a}\,dr &&: a>r\\
&\int_0^a-\frac{4\pi G r^2\rho(r)}{a}\,dr+\int_a^R-4\pi G \rho(r) r\,dr &&: a\leq r
\end{aligned}
\right.
\\
=\left\{\begin{aligned}
&-\frac{GM(R)}{a} &&: a>r\\
&-\frac{GM(a)}{a}+\int_a^R-4\pi G \rho(r) r\,dr &&: a\leq r
\end{aligned}
\right..$$
In essence, your expression for potential was missing the term $\int_a^R-4\pi G \rho(r) r\,dr$, which caused the errors you saw.
Finally, answering your question
We then compute the gradient, $\nabla\Phi(a)$, using the $a\leq R$ case for $\Phi$ listed above. In Mathematica this easily becomes:
-D[Integrate[-((4*G*Pi*r^2*\[Rho][r])/a), {r, 0, a}] + 
  Integrate[-4*G*Pi*r*\[Rho][r], 
       {r, a, R}], a]//TeXForm

yielding
$$a_g(a)=-\nabla\Phi(a)=-\int_0^a \frac{4 \pi  G r^2 \rho (r)}{a^2} \, dr=-\frac{GM(a)}{a^2},$$
exactly as you wanted.
Lane-Emden in a nutshell
In essence, Lane-Emden is just the combination of Poisson's equation $\Delta\Phi=4\pi G\rho$ (which is universally valid), Newton's law $\rho\nabla\Phi=\nabla P$ (which is valid in a gravitating fluid), the polytropic constitutive relation $P=C\rho^k$ (which is an approximation), and spherical symmetry. The polytropic assumption allows you to eliminate $P$, resulting in a differential equation solely in $\rho$.
(Minor note: Lane-Emden is commonly expressed in a dimensionless parameter $\theta$, rather than $\rho$ itself.)
