Why is the quadrupole moment of a spherical object equal to zero (when looking at the formula in cartesian coordinates)? Is there a simple way to understand why
$$\int \rho(x,y,z) (2z^2-x^2-y^2) dxdydz$$
is equal to zero if the density has spherical symmetry?
 A: Spherical symmetry makes $x$, $y$, and $z$ “equivalent directions” and $1+1-2=0$. (An integrand of $x^2$, $y^2$, or $z^2$ should produce the same integral when $\rho$ is spherically symmetric.)
A: If it's spherically symmetric then you can interchange or invert x,y,z and have the same density, e.g. $\rho(x,y,z) = \rho(-y,z,x)$
This causes the second part of the integral $(2z^2 - x^2 - y^2)$ to average out zero for the same density.
The integral has contributions from 8 non-overlapping octants plus the boundaries between. Because $\rho$ is spherically symmetric the integral receives equal contributions from all 8 octants. W.L.O.G. consider the octant where all coordinates are positive.
Because $\rho$ is spherically symmetric we can interchange any coordinates and get the same value e.g. $\rho(x,y,z) = \rho(z,x,y)$ (this corresponds to some mirroring and rotation of the sphere). Now split the octant into 6 parts depending on coordinate order (e.g. $x < y < z$). Take some point in one of these parts and consider the contributions to the integral from the corresponding points in all 6 parts (i.e. all 6 re-orderings of their coordinates). The density function is the same at all 6 points because of the spherical symmetry. Because of the second term, the contribution to the integral from these 6 points cancels itself out:
$$
(2z^2 - x^2 - y^2) +\\(2x^2 - y^2 - z^2) +\\(2y^2 - z^2 - x^2)+\\
(2z^2 - y^2 - x^2) +\\(2x^2 - z^2 - y^2) +\\(2y^2 - x^2 - z^2)\\
= (2-1-1+2-1-1)z^2\\
+ (-1-1+2-1-1+2)y^2\\
+ (-1+2-1-1-1+2)x^2\\
= 0z^2 + 0y^2 + 0x^2\\
= 0
$$
This holds for all points inside such a segment of an octant, therefore when considered together with their corresponding points in other segments, none of the points any segment contribute anything to the integral, but that's every point, therefore there is no contribution to the integral at all and the integral is zero.
As for why the boundaries between segments and between octants (where two coordinates are equal or one of them is zero) can't add up to something and make the integral nonzero, I handwave them away by pointing out that 0% of the points lie on the boundaries and therefore they contribute 0% to the integral. You did ask for simple, rather than rigorous.
A: Ghoster has already given the shortest possible answer, demonstrating how the result can be seen without any calculation within less than a millisecond.
Hopefully not ending up with the longest possible answer, let me discuss some underlying group theoretic aspects, relevant in many physical problems where invariance under rotations plays a crucial role.
Consider the second rank tensor $$T_{ij} = \int\limits_{\mathbb{R}^3} \! d^3 \!x \, \rho(|\vec{x}|)\,  x_i x_j,$$ where $\rho(|R \vec{x}|) = \rho(|\vec{x}|)$ for arbitrary rotations $R \in \rm SO(3)$, i.e. $R^T=R^{-1}$ and $\det R =1$. $T_{ij}$ is symmetric under an interchange of the indices $ i \leftrightarrow j$ and also invariant under rotations, as $$T^\prime_{ij}= R_{ik} R_{j \ell} T_{k\ell}= \int\limits_{\mathbb{R}^3}\! d^3 \! x \, \rho(|\vec{x}|)(R_{ik}x_k) (R_{j  \ell}x_{\ell}) =  \int\limits_{\mathbb{R}^3}\! d^3 \! y \, \rho(|R^{-1} \vec{y}|) \, y_i y_j = \int\limits_{\mathbb{R}^3} \! d^3 y \, \rho(|\vec{y}|) \, y_i y_j =T_{ij},$$ where $\det R =1$ was used for the Jacobian of the change of variables $\vec{y}= R \vec{x}$ in the integral. As a consequence of Schur's lemma, a symmetric invariant tensor of rank two must be proportional to the Kronecker delta, $T_{ij}= \lambda \delta_{ij}$. This can easily be seen by rewriting the relation $R_{ik} R_{j \ell} T_{k \ell} = T_{ij}$ in matrix form, $R \, T R^T=T \Leftrightarrow [R, T]= 0$. As $T$ commutes with all elements $ R \in \rm SO(3)$, it must be a multiple of the identity matrix, $T= \lambda \mathbb{1}_3$.
Contracting both indices in $$ T_{ij}= \lambda \delta_{ij} = \int\limits_{\mathbb{R}^3} \! d^3 x \, \rho(|\vec{x}|) \, x_i x_j,$$ we find $$\lambda = \frac{1}{3} \int\limits_{\mathbb{R}^3} \! d^3 x \, \rho(|\vec{x}|) \, \vec{x}^2 = \frac{4 \pi}{3} \int\limits_0^\infty \! dr \, r^4 \rho(r), $$ yielding the final result $$ T_{ij} = \int\limits_{\mathbb{R}^3} \! d^3 x \, \rho(|\vec{x}|) \, x_i x_j= \frac{4 \pi}{3} \int\limits_0^\infty \! dr \ r^4 \rho(r) \, \delta_{ij},$$ thus reducing the determination of $6$ three-dimensional integrals to a single one-dimensional integration. (The integral of OP is just $2 T_{33}-T_{22}-T_{11}= \lambda (2-1-1)=0$.)
Also higher rank tensor integrals with a rotation invariant measure can be treated in a similar fashion. As an example, the totally symmetric invariant tensor $$S_{ijk\ell}= \int\limits_{\mathbb{R}^3} \! d^3 x \, \rho(|\vec{x}|) \, x_i x_j x_k x_\ell $$ takes the form $S_{ijk\ell}= c (\delta_{ij} \delta_{k\ell}+\delta_{ik}\delta_{j\ell}+\delta_{i\ell}\delta_{jk})$ with $c= \frac{4 \pi}{15}\int\limits_0^\infty \! dr \, r^6  \rho(r)$.
