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.
First, $\vec{r}^\prime$ is a vector that goes from the origin to the source of charge. If the source is a volumetric distribution, one must sum all contributions of charge, that's why one integrates over all the volume, say $\mathcal{V}$; the (correct) expression for the potential should be
$$V( \vec{r}) = \frac{1}{4 \pi \epsilon _{0}} \int_\mathcal{V} \frac{\rho (\vec{r}^\prime)}{ℛ} d\mathcal{V}^\prime$$
so that all dependence of $V$ remains on $\vec{r}$. Then, $r^\prime$ is just the magnitude $|\vec{r}^\prime|$, being the distance from the origin to the source of charge.
Second, usually, the series expansion of a function $f(x)$ about some point $x_0$ is useful because if you want to know the value of $f$ near $x_0$, you may just take some few terms of the expansion; it is as seeing the plot of $f$ with a magnifying glass. You should remember this from your first calculus courses, it is done a lot in physics. Here the expansion about $\epsilon=0$ will be useful since $\epsilon\to0$ implies $r\to\infty$ (just really big, if you will). The (correct) expression
$$V(\vec{r}) = \frac{1}{4 \pi \epsilon _{0}} \sum ^{\infty}_{n=0}\frac{1}{r^{n+1}} \int(r')^n\,P_{n}(\cos \theta^\prime)\,\rho( \vec{r}')\,d\mathcal{V}'$$
is just another way of writing the series expansion in terms of $r$, $r^\prime$ and $\theta^\prime$, where $P_n$ are the Legendre polynomials (Griffiths defines them there, ain't he?). This expression is useful, as it means, explicitly, that
$$V(\vec{r})=\frac{1}{4\pi\epsilon_0}\left[\frac{1}{r}\int\rho(\vec{r}')\,d\mathcal{V}'+\frac{1}{r^2}\int{r'}\cos\theta'\,\rho(\vec{r}')\,d\mathcal{V}'+\frac{1}{r^3}\left(\cdots\right)+\ldots\right]$$
so that if you want to evaluate the potential for points far from the source (big $r$), then you may just neglect higher order terms in $r$ and just take the $1/r$ (monopole) term; and so on if you're considering a better approximation, you may take the $1/r^2$ (dipole) term, etc... That's the real usefulness of the series expansion; in a lot of situations evaluating $V( \vec{r}) = \frac{1}{4 \pi \epsilon _{0}} \int \frac{\rho (\vec{r}')}{ℛ} d\mathcal{V}'$ will get really ugly, and then, mostly, is when the multipole approximation will be useful.
Best Answer
The multipole expansion in its most general form reads $$ V(\mathbf r)=\sum_{l=0}^\infty \sum_{m=-l}^lQ_{lm}\frac{Y_{lm}(\theta,\phi)}{r^{l+1}}, $$ where the $Q_{lm}$ are the multipole moments of the system and the $Y_{lm}$ are spherical harmonics, and in this form it is applicable to bodies of any shape, and located at any point in space, so long as the body is contained in a finite spherical region of radius $R$ and the test point $\mathbf r$ is located outside this spherical region.
The Wikipedia article you linked to is concerned with the first few terms in this expansion, which can be rewritten as
$$ V(\mathbf r)=-\frac{GM}{r}+\frac{\mathbf d\cdot\mathbf r}{r^3}+\frac{p_2(x,y,z)}{r^5}, $$ where $p_2$ is a traceless homogeneous quadratic polynomial, and the dipole moment of the system is defined as $$ \mathbf d=G\int\mathbf r\:\mathrm dm(\mathbf r). \tag 1 $$
These first three terms go down with the distance to the origin $r$ as $1/r$, $1/r^2$ and $1/r^3$ respectively, which means that the relative importance of the second, dipole term goes as $1/r$ with respect to the first, monopole term.
This changes, however, if you set up your origin at the centre of mass of the distribution, since then the dipole moment $(1)$ vanishes. If you do this, then you can completely ignore the dipole term, and the subleading term of the expansion becomes the quadrupole term, which goes as $1/r^2$ of the monopole term, which makes the truncated series much more accurate.
That's about it, really.
In addition, though, you should note that the centre of gravity as you've linked to cannot be defined for a body in isolation: it is a property of a mass distribution $\mathrm dm(\mathbf r)$ and an externally-produced gravitational field $\mathbf g(\mathbf r)$ acting on the distribution. In these conditions you can find the total force $$ \mathbf F=\int\mathbf g(\mathbf r)\mathrm dm(\mathbf r) $$ and the total torque with respect to an arbitrary origin $\mathbf r_0$, $$ \boldsymbol\tau=\int(\mathbf r-\mathbf r_0)\times\mathbf g(\mathbf r)\mathrm dm(\mathbf r), $$ and then you can define the centre of gravity as the $\mathbf r_0$ at which $\boldsymbol \tau$ vanishes. If you don't have the external field $\mathbf g(\mathbf r)$, though, this doesn't make any sense: the centre of gravity is the point through which a given field configuration can be thought of as acting, and it depends on the field configuration.