This is a delightfully interesting question, and @DavePhD has put up an excellent answer. I will try and contribute my own perspective to it, although more or less in line with what has already been said; I have also included some illustrations to make concepts a bit more tangible, and perhaps, palatable
Preliminaries
How do we determine if something is spherical? Well, we do so by "looking at", measuring its dimensions etc. For large macroscopic objects, that usually means photons of light scatter off the object, our ruler etc. and are captured by our eyes. The signals generated are converted into visual information by our brains. Essentially, all processes involved are quantum mechanical, however, for macroscopic objects their unintuitive, quantum weirdness does not manifest itself.
Similar, without going into details, to determine the "spherical-ness" of an atom we need to probe some physical observable, with which we associate an operator $\hat{\mathcal{O}}$. For a given quantum state $|\psi\rangle$, the expectation value for this observable, i.e. something we expect to get would be given by $\langle \psi|\hat{\mathcal{O}} |\psi\rangle$ (depending on your level of education, you may need to look this up. You can read some details here). Anyway, we can also talk about the probability density of finding an electron in some region of space, and that would be given by
$\langle \psi|\psi\rangle \equiv |\psi|^2 $
So essentially, one never really sees the "wavefunction".
Also, one needs to point out is that space is isotropic. Basically, this means there is nothing special about any particular direction in three-dimensional space. Also, even from classical electrostatics we know that Coulombic interactions behave the same way: if you have a negative charged place at some point in space, it will attract positive charges all around it in the same fashion regardless of whether they are placed to its right or two its left.
Now, we know electrons in atoms occupy "orbitals" which are obtained by solving Schrödinger's equation. It is complicated to solve for many-electron systems, but we can do it for hydrogen. Writing down the equation, we see,
$\left( -\frac{\hbar^2}{2\mu} \nabla^2 - \frac{e^2}{4\pi \epsilon_o}\right)\psi (r,\theta, \phi) = E \psi (r,\theta, \phi) $
Here, $\mu$ is the reduced mass of the system, this allows us to treat the problem in centre of mass coordinates (remember Coulomb forces are central). $\nabla^2$ is the laplacian operator, which basically contains the second-derivatives along x,y,z in cartesian space. It has a more complicated form in spherical coordinates which is what we are using for our wavefunction $\psi$, but basically, the important fact is that this first term is kinetic energy, and then the second term is the coulombic potential due to the interaction of electron and nucleus.
This is a separable, partial differential equation which can be solved in terms of special functions. The normalised position wave functions, given in spherical coordinates are
$\psi_{nlm} (r,\theta, \phi) = \sqrt{\left(\frac{2}{na_0^*}\right)^3 \left(\frac{(n-l-1)!}{(2n+l)!}\right)} e^{-\rho/2}\rho^l L^{2l+1}_{n-l-1} (\rho) Y^m_l( \theta, \phi) $
Where, $\rho = \frac{2r}{na_{0}^{*}}$,
$a_{0}^{*} = \frac{4 \pi \epsilon _{0}\hbar ^{2}}{\mu e^{2}}$ (reduced Bohr radius)
$L$ and $Y$ denote the generalised Laguerre Polynomial and Spherical Harmonic respectively, and n,l,m are the quantum numbers you must be familiar with.
Anyway, we are interested only in the functional forms, so we can ignore all the constants out in front. We also note that the wave function is a product of two-functions, radial and angular.
$\psi_{nlm} (r,\theta, \phi) = R(r)Y^m_l(\theta, \phi)$
Now, essentially, spherical symmetry demands that there will be no angular dependance, and any properties measured or probability of finding electrons etc. only depends on distance from the centre r. This means whenever you measure $\langle \psi|\hat{\mathcal{O}} |\psi\rangle$, it will again only depend on $r$ and not on the angular variables.
We only need to focus our attention on the Spherical Harmonics, as the contain the angular information.
For the purpose of this discussion, we can crudely assume that even in many-electron atoms, the spatial wave functions are similar to hydrogenic wave functions (which is only a one electron system). This is because it is only possible to obtain analytic, functional solutions for the hydrogen atom.
Main Argument
For say an $s$-orbital like l = 0, m = 0. It is trivial to see that the wave function, and consequently $|psi|^2$ is spherically symmetric. i.e. depends only on $r$ because $Y^0_0 = \frac{1}{2\sqrt{ \pi}}$.
Now, let us look at the $p$-orbitals where $l = 1$ and $m$ takes values $-1,0,1$
$$ Y^0_1 = \frac{1}{2} \sqrt{\frac{3}{\pi}} cos \left( \theta \right)$$
$$ Y^{\pm 1}_1 = \mp \frac{1}{2} \sqrt{\frac{3}{\pi}} sin \left( \theta \right)e^{\pm i \phi}$$
Now, clearly $\sum\limits_{j} |Y^j_1|^2 = \frac{3}{4\pi}$ is independent of $\theta$ and $\phi$
The same analysis can be repeated for $d,f$ -orbitals as well. A table of spherical harmonics can be found here
So essentially atoms with electrons in valence $s$-orbitals and/or completely filled or half-filled $p,d,f$ shells would be spherically symmetric.
Now, what if we have something like $2p^1$ type situation. A single $p$-orbital is not spherical. But given that space is isotropic, we have no external fields etc., all the $p$-orbitals corresponding to different $m$ values are essentially the same (degenerate). Nature can't tell the difference, and neither can you when you perform an experiment to measure the "Spherical-ness" of an electron.
This does not mean we are naively taking a superposition of the $p$-orbitals. If you see proof by pictures section, this is clearly not true
Proof by Pictures
Note: The spherical harmonics are complex valued functions, so they are hard to plot, the individual $p$ orbitals conventionally plotted are generated by taking linear combinations of them which are still solutions to the Schrödinger's equation.
Given below is $p_z$ orbital, $p_x$ and $p_y$ look the same just are oriented along the $x$ and $y$ axis.

A naive sum of $p_x$, $p_y$ and $p_z$ orbitals.

Plot of $\sum\limits_{j} |Y^j_1|^2 $. Remember that observable quantities are related to $|\psi|^2$.
