Short answer
The results come from:
- the assumption of homogeneity (no preferential point in space) and isotropy (no preferential orientation) of space,
- a reasonable choice of parameters to describe such a space (spherical coordinates in 3D space, cylindrical in 2D space),
- an application of divergence theorem to proof the result does not depend on the domain of integration.
Let's answer the points in your question.
Poisson equation, active/passive points
But before focusing on these points, let's inspect Poisson equation with Dirac's delta force $\delta(\mathbf{x} - \mathbf{x}_0)$
\begin{equation}
-\nabla^2 G(\mathbf{x}, \mathbf{x_0}) = \delta(\mathbf{x} - \mathbf{x}_0) \ ,
\end{equation}
and give a meaning to $\mathbf{x}$ and $\mathbf{x}_0$ points in space:
- point $\mathbf{r}_0$ can be interpreted as the point where we can find the "singularity" of the Dirac delta, and it can be interpreted as a parameter of the problem; since it's the location of the forcing, let's call it active point;
- point $\mathbf{x}$ can be interpreted as the independent variable of the problem, i.e. the set of point in space $\mathbf{x} \in V$ where we want to evaluate the function $G$; let's call it passive point.
Parametrization
Following the assumption of space homogeneity and isotropy, it's reasonable to look for a solution with spherical symmetry around the active point $\mathbf{x}_0$. Let's thus:
- introduce a set of spherical coordinates with origin in $\mathbf{x}_0$, so that
- each point in space can be written as $\mathbf{x} = \mathbf{x}_0 + \mathbf{r}$, being $\mathbf{r}$ the vector connecting $\mathbf{x}$ and $\mathbf{x_0}$, the origin of the spherical coordinate system,
- the radial coordinate of a point $\mathbf{x}$ reads $r = |\mathbf{r}| = |\mathbf{x} - \mathbf{x}_0|$
- assume that the solution of the problem, the Green's function, only depends on $r$ (and thus on $|\mathbf{r}| = |\mathbf{x} - \mathbf{x}_0|$), and not on the angles $\phi$, $\theta$.
I. Integration volume (I). Performing an integration over a spheres and its surface of radius $R$, the author identifies all the points of the sphere $\mathbf{x}$ with their distance from the origin: 2-dimensional surface in 3d-space only needs an equation (in simple cases, like a sphere) to be determined; thus, instead of writing $\mathbf{x} \in S_R$, he identifies this points with the distance from the origin, i.e. the radial coordinate $r=R$.
II. Integration volume (II). The result doesn't depend on the choice of the domain of integration, but only on the assumptions of homogeneity and isotropy. You can prove it using divergence theorem on a "shell" volume $\Omega$ that doesn't contain point $\mathbf{r}_0$, so that Poisson equation in $\Omega$ is homogeneous
\begin{equation}
-\nabla^2 \mathbf{G}(\mathbf{x},\mathbf{x}_0) = 0 \ .
\end{equation}
\begin{equation}
\begin{aligned}
0 & = \int_{\Omega} \nabla^2 \mathbf{G} = \\
& = \oint_{\partial \Omega} \nabla G \cdot \mathbf{\hat{n}} = 0 \\
& = \oint_{S_R} \nabla G \cdot \mathbf{\hat{n}} + \oint_{\tilde{S}} \nabla G \cdot \mathbf{\hat{n}}^- \ ,
\end{aligned}\end{equation}
where surfaces $S_R$ and the arbitrary surface $\tilde{S}$ are the outer and the inner surface of the "shell" volume $\Omega$. Here, surface $\tilde{S}$ is supposed to be the inner surface of the shell, and thus with reversed unit normal $\mathbf{\hat{n}}^- = \mathbf{\hat{n}}$.
Thus, it's easy to realize that
\begin{equation}
\oint_{S_R} \nabla G \cdot \mathbf{\hat{n}} = \oint_{\tilde{S}} \nabla G \cdot \mathbf{\hat{n}} \ .
\end{equation}
A spherical domain is usually chosen, in order to get integral to evaluate as simple as possible.
III. Space and parametrization, $\mathbf{r}$, $\mathbf{x}$, $r$, $R$. Once you got the result for all the points in space lying on the sphere of radius $R$ centered in the Dirac's delta location, and determined by the radial coordinates $r=R$,
\begin{equation}
G(R) = - \dfrac{1}{4\pi R} \ ,
\end{equation}
this solution holds for all the points in space, and thus for every value of the radial coordinate $r$,
\begin{equation}
G = -\dfrac{1}{4\pi r} = -\dfrac{1}{4 \pi |\mathbf{r}|} = -\dfrac{1}{4\pi |\mathbf{x} - \mathbf{x}_0|} \ .
\end{equation}
Last remark on notation
To be honest and very precise, when you change the independent variables, you should also change the name or the letter you use to indicate the functions, since mathematically these functions different relations (roughly speaking, a function is nothing more than a relation between variables) between independent and dependent variables, i.e.
\begin{equation}
G(r) = -\dfrac{1}{4 \pi r} = G(|\mathbf{r}|) = -\dfrac{1}{4 \pi |\mathbf{r}|} = G(|\mathbf{x}-\mathbf{x}_0|) = -\dfrac{1}{4 \pi |\mathbf{x}-\mathbf{x}_0|} = \tilde{G}(\mathbf{x}, \mathbf{x}_0) \ ,
\end{equation}
being the function $\tilde{G}(\mathbf{a}, \mathbf{b})$ the composition of the function $N(\mathbf{a}, \mathbf{b}) = |\mathbf{a} - \mathbf{b}|$ and $G(r)$,
\begin{equation}
\tilde{G}(\mathbf{a}, \mathbf{b}) = G(|\mathbf{a} - \mathbf{b}|) = G(N(\mathbf{a}, \mathbf{b})) \ .
\end{equation}