Solution 1:

The first two observations can be explained by the fact that if $r$ is the root of $a_nz^n+\cdots +a_0$ then $r^{-1}$ is the root of $a_0z^n+\cdots +a_n$. Since the joint probability density of the coefficients is symmetric under $a_k\mapsto a_{n-k}$ transformation the density of the roots should be symmetric under $r\mapsto r^{-1}$. It should be fairly straightforward to find asymptotic behavior of the density as $z\to0$ or $z\to\infty$. In the first case, $r\to a_0/a_1$ and the second case is symmetric under inversion.

The distribution around $|z|=1$ is less obvious, so I need to think more to explain your third observation.