Where am I going wrong in deriving Maxwell's first equation in differential form? By denoting the source coordinates with prime, I get flux through a closed surface:
$$\Phi= 
\displaystyle\oint_{A}  \mathbf{E}(x,y,z) \cdot \mathbf{\hat{n}}\ dA 
=q (x',y',z')$$
And now using the formal definition of divergence:
$$\displaystyle \nabla \cdot \mathbf{E}(x,y,z)=\dfrac{d\Phi}{dV} = \dfrac{d^3\ q (x',y',z')}{dx\ dy\ dz}=0 \tag1$$
In order to correctly get $\nabla \cdot \mathbf{E}(x,y,z)=\rho (x,y,z)$, we should be having:
$$\nabla \cdot \mathbf{E}(x,y,z)=\dfrac{d^3\ q(x,y,z)}{dx\ dy\ dz} \tag2$$
instead of equation $(1)$. Where am I going wrong?

EDIT: A simpler way of asking the same question:
By divergence theorem:
\begin{align}
\iiint \nabla \cdot \mathbf{E}\ dV &= \unicode{x222F} \mathbf{E} \cdot \hat{\mathbf{n}}\ dS \\
&=q\ (x',y',z') \\
&=\iiint \rho (x',y',z')\ dV'
\end{align}
Now how is the cancellation of triple integral in LHS and RHS justifiable when LHS has integral with respect to $V$ while RHS has integral with respect to $V'$?

 A: Just because the two integrals over the same volume are equal, that does not imply that the integrands are equal.  However, if the integral relation
$$ \int _V \nabla \cdot \mathbf{E} \ dV = \int _V \rho \ dV$$
holds for any volume, then it can be applied to an infinitesimally small volume, $V_{i}$, for which $\rho$ is a constant.  Then you get the relation
$$ \int _{V_i} \nabla \cdot \mathbf{E} - \rho \ dV = 0$$
This can be applied everywhere, pointwise, and yields the differential form.
A: It looks like you are using the convention where unprimed coordinates represents The point of space you are "interested in" and primed coordinates refer to the coordinates of the charge density. So, for example, Coulomb's law becomes
$$\mathbf E(\mathbf r)=\iiint\frac{\rho(\mathbf r')(\mathbf r-\mathbf r')}{|\mathbf r-\mathbf r'|^3}\ \text dV'$$
However, for what you are trying to do the distinction between primed and unprimed coordinates is irrelevant. This is because for Gauss's law in integral form that you start with we don't need to distinguish between these two spatial coordinates. You're defining a surface and looking at the field on that surface and the charge contained within that surface. Therefore you only have one set of spatial coordinates to worry about, not two. 
It might be better for you to consider a single set of coordinates for something like Coulomb's law anyway. We specify a fixed point in space $\mathbf r_0= x_0\hat x+y_0\hat y+z_0\hat z$ and then integrate over the spatial coordinates containing our charge distribution:
$$\mathbf E(\mathbf r_0)=\iiint\frac{\rho(\mathbf r)(\mathbf r_0-\mathbf r)}{|\mathbf r_0-\mathbf r|^3}\ \text dV$$
We are really only using one set of coordinates here along with a specified member of these coordinates.
So, putting together your whole derivation here, we start with Gauss's law in integral form:
$$\displaystyle\oint\mathbf E(\mathbf r)\cdot\text d\mathbf A=q_{enc}$$
where we are integrating over an arbitrary surface on the left and the right is how much charge in contained within this surface. Note that the enclosed charge is not a function, it is a single value based on the surface we have chosen. Also notice that the infinitesimal area vector $\text d\mathbf A$ and the vector $\mathbf r$ refer to the same spatial coordinates.
Then we use the divergence theorem on the left and the definition of charge density on the right:
$$\iiint\nabla\cdot\mathbf E(\mathbf r)\ \text dV=\iiint\rho(\mathbf r)\ \text dV$$
Notice that these two integrals are over the same volume, which is the volume contained by the surface we chose earlier. As has been mentioned in other answers, since the surface we chose was arbitrary we can say these integrands are equivalent.
A: Note that $q$, in general, is not a function of space coordinates. The density $\rho$ is. So the RHS of the Maxwell equation really is
$$
\int \rho(x, y, z) dV
$$
