For question 2: ("Why does a single charge away from the origin have a dipole term?")
Let's say you have a charge of +3 at point (5,6,7). Using the superposition principle, you can imagine that this is the superposition of two charge distributions
Charge distribution A: A charge of +3 at point (0,0,0)
Charge distribution B: A charge of -3 at point (0,0,0) and a charge of +3 at (5,6,7).
Obviously, when you add these together, you get the real charge distribution:
$$
(\text{real charge distribution}) = (\text{charge distribution A}) + (\text{charge distribution B}).
$$
By the superposition principle:
$$
(\text{Real }\mathbf E\text{ field}) = (\mathbf E\text{ field of charge distribution A}) + (\mathbf E\text{ field of charge distribution B}).
$$
And, since the multipole expansion also obeys the superposition principle:
\begin{align}
(\text{real monopole term}) & = (\text{monopole term of distribution A}) + (\text{monopole term of distribution B}),\\
(\text{real dipole term}) & = (\text{dipole term of distribution A}) + (\text{dipole term of distribution B}),\\
(\text{real quadrupole term}) & = (\text{quadrupole term of distribution A}) + (\text{quadrupole term of distribution B}),
\end{align}
and so on.
The field of charge distribution A is a pure monopole field, while the field of charge distribution B has no monopole term, only dipole, quadrupole, etc. Therefore,
\begin{align}
(\text{real monopole term}) & = (\text{monopole term of distribution A}), \\
(\text{real dipole term}) & = (\text{dipole term of distribution B}),\\
(\text{real quadrupole term}) & = (\text{quadrupole term of distribution B}),
\end{align}
and so on.
Even though it's unintuitive that the real charge distribution has a dipole component, it is not at all surprising that charge distribution B has a dipole component: It is two equal and opposite separated charges! And charge distribution B is exactly what you get after subtracting off the monopole component to look at the subleading terms of the expansion.
I would say the quickest route would be to use the multipole expansion for the potentials of the different charges and then add those up. Thus you'd use
$$
\frac{1}{|\mathbf{r}-\mathbf{r}'|}=\sum_{l=0}^\infty\sum_{m=-l}^l\frac{ 4\pi}{2l+1} \frac{r'^l}{r^{l+1}}Y_{lm}^\ast(\theta',\phi')Y_{lm}(\theta,\phi),
$$
valid for $r>r'$ (see Jackson, 2nd ed., p.111, eq 3.70). To implement this, take $\mathbf{r}'=(0,0,\pm d)$, so that $\theta'=0,\pi$ and $\phi'$ is irrelevant. Then $Y_{lm}^\ast(\theta',\phi')$ is zero unless $m=0$, so you have
$$
Y_{lm}^\ast(\theta',\phi')=\delta_{m,0}N_{lm}P_l(\cos\theta)=\delta_{m,0}\sqrt{\frac{2l+1}{4\pi}}P_l(\pm1)=\delta_{m,0}\sqrt{\frac{2l+1}{4\pi}}(\pm1)^l.
$$
Thus, in this case,
$$
\frac{1}{|\mathbf{r}-\mathbf{r}'|}=\sum_{l=0}^\infty (\pm d)^l \sqrt{\frac{ 4\pi}{2l+1}} \frac{Y_{l0}(\theta,\phi)}{r^{l+1}}=\sum_{l=0}^\infty (\pm d)^l \frac{P_{l}(\cos\theta)}{r^{l+1}}.
$$
If you put it all together, then, the $l=0$ components of the charges at $(0,0,\pm d)$ will cancel with the one at the origin (all the other multipole terms of which vanish, of course), and the $l=1$ terms will kill each other due to that sign. The full potential will therefore be
$$
\Phi=q\sum_{l=2}^\infty \left[1+(-1)^l \right] d^l \frac{P_{l}(\cos\theta)}{r^{l+1}}=\sum_{k=0}^\infty 2q d^{2+2k} \frac{P_{2+2k}(\cos\theta)}{r^{3+2k}}.
$$
Now, as to how to extract the multipole moments, it depends on what normalization you use. I would use a normalization as
$$
\Phi=\sum_{l=0}^\infty\sum_{m=-l}^l q_{lm}\frac{Y_{lm}(\theta,\phi)}{r^{l+1}},
$$
but it depends on what the multipolar fields are defined like in your book. Since the nonzero-$m$ terms have already vanished, we must have $q_{2,\pm1}=q_{2,\pm2}=0$ as the problem's symmetry demands (as I indicated already in the comments). The $m=0$ component is then simply
$$q_{2,0}=2q d^2,$$
but, again, this depends on the normalization you're using.
Best Answer
The full form of the expression for $Q_0$ is $$ Q_0=\frac{1}{e}\int_0^{2\pi}\mathrm d\phi\int_0^r r^4\mathrm dr\int_0^\pi (3\cos^2(\theta)-1)\sin(\theta)\:\rho(r,\theta,\phi)\:\mathrm d\theta. $$ If you assume that $\rho(r,\theta,\phi)$ has no angular dependence, i.e. that your nucleus is spherical, then indeed the inner integral reduces to $$ \int_0^\pi (3\cos^2(\theta)-1)\sin(\theta) \:\mathrm d\theta=0, $$ giving you $Q=0$. If you do not makes such an assumption, you cannot get rid of the $\theta$ dependence of $\rho$ in that inner integral, which can make it nonzero.
In other words, yes: spherical nuclei are spherical. Non-spherical nuclei are not.