PRIMER:
Note that in regions that are absent of charge, $\nabla \cdot E=0$. Let $V$ be a region, bounded by the surface $S$, that contains a point charge. Furthermore, let $V'$ be a region bounded by the surface $S'$, contained in $V$.
If the point charge is in $V'$, then it is outside of the region $V\setminus V'$ and from Gauss's law and the Divergence Theorem we have
$$\begin{align}
\int_{V\setminus V'}0\,dV&=\int_{V\setminus V'}\nabla\cdot\vec E\,dV\\\\
&=\oint_{S} \vec E\cdot \hat n\,dS-\oint_{S'} \vec E\cdot \hat n\,dS
\end{align}$$
from which we conclude
$$\oint_{S} \vec E\cdot \hat n\,dS=\oint_{S'} \vec E\cdot \hat n\,dS$$
That is, the value of the surface integral is independent of the surface over which the integration is taken.
If the point charge is in $V\setminus V'$, then $\nabla \cdot E=0$ in $V'$ and $\oint_{S'}\vec E\cdot \hat n\,dS=0$.
We now will show explicitly that if $S'$ is any sphere that contains a point charge, that $\oint_{S'}\vec E\cdot \hat n\,dS=q/\varepsilon_0$ (We could have even chosen a sphere for which the point charge is at the origin.). To that end we proceed.
The electric field $\vec E(\vec r')$ due to a point charge at $\vec r$ at points on the surface of a sphere, centered at $0$ and with radius $a$ is
$$\vec E(\vec r')=\frac{q(a\hat r'-\vec r)}{4\pi\varepsilon_0|\vec r-a\hat r'|^3}$$
So, we wish to compute the integral
$$\oint_{|\vec r'|=a} \vec E(\vec r')\cdot \hat n'\,dS'=q\int_0^{2\pi}\int_0^\pi \left(\frac{a\hat r'-\vec r}{4\pi \epsilon_0|\vec r-a\hat r'|^3}\right)\cdot\hat r'\,a^2\,\sin(\theta')\,d\theta'\,d\phi'$$
Due to the spherical symmetry, without loss of generality, we align the polar axis so that the point charge lies at $\vec r=r\hat z$. Then,
$$\begin{align}
\oint_{|\vec r'|=a} \vec E(\vec r')\cdot \hat n'\,dS'&=\frac{qa^2}{2\varepsilon_0}\int_0^\pi \frac{a-r\cos(\theta')}{(r^2+a^2-2ar\cos(\theta'))^{3/2}}\sin(\theta')\,d\theta'\\\\
&=\frac{qa^2}{2\varepsilon_0}\int_{-1}^1 \frac{a-rx}{(r^2+a^2-2arx)^{3/2}}\,dx\tag1
\end{align}$$
The integral on the right-hand side of $(1)$ is
$$\begin{align}\int_{-1}^1 \frac{a-rx}{(r^2+a^2-2arx)^{3/2}}\,dx&=-\left.\left(\frac{r-ax}{a^2(r^2+a^2-2arx)^{1/2}}\right)\right|_{-1}^1\\\\
&=\frac{r+a}{a^2(r+a)}-\frac{r-a}{a^2|r-a|}\\\\
&=\begin{cases}
0&, r>a\\\\
\frac 2{a^2}&,r<a\tag2
\end{cases}
\end{align}$$
Using $(2)$ in $(1)$ we find the coveted result
$$\oint_{|\vec r'|=a} \vec E(\vec r')\cdot \hat n'\,dS'=\begin{cases}\frac{q}{\varepsilon_0}&,r<a\\\\
0&,r>a\end{cases}$$
as was to be shown!