In my discussion here I am assuming your $\theta$ is effectively a longitude and $\phi$ is effectively a latitude. Perhaps more typical spherical co-ordinates use an angle down from the north pole rather than up from the equator and swap the roles of the two symbols from that - but it's no problem to deal with it either way, so I'll stick with what your notation appears to be.
Note that the distribution of radius is not of interest here, only the angles, so we could project everything onto a unit sphere without changing the angles. This is quite useful in the normal case.
With a spherically symmetric distribution like the three dimensional standard normal, the appearance of the distribution of inclination has to do with the fact that there's a lot more area on the surface of a sphere near the equator than near the poles.
If you follow through the mathematics (or write a geometric argument in terms of elements of probability similar to the earlier 2D question), you can get that the inclination should have a density proportional to $\cos(\phi)$. Here's a geometrical argument that should motivate it in the "elements of probability" terms:
Since the radius at the equator is 1 and the radius at latitude $\phi$ is $\cos(\phi)$, the circumference at latitude $\phi$ is proportional to $\cos(\phi)$, and so the density at $\phi$ is proportional to $\cos(\phi)$.
Uniform case: With the 3D-uniform normalized to constant radius, you don't have uniformity of density on the sphere for the same reason that we didn't in the 2D case - when you project onto the sphere, there's a lot more "density" on the sphere near the angles where the corners are than where the sides are (with parts near the middle of the edges being in between) -- because there's more of the volume of the cube for angles close to the corners than for angles near the middle of the faces.
We can see this by generating many random values uniformly in the cube and projecting them onto the sphere. Since there's more volume near the corners than near the faces of the cube, there's greater density looking "inward" from the corners than the faces. If we plot the height (recall this is a projected-z-value, $z^* = z/r$, where $r=\sqrt{x^2+y^2+z^2}$) above the equator against the longitude, we get the top plot below:
That height corresponds to the vertical side of the right triangle in the previous diagram; that height is the $\sin$ of $\phi$ ($z^*=\sin(\phi)$). To convert that to the latitude ($\phi$), we would take the arcsin of that projected vertical height, which is what we see in the lower plot. This "stretches" things more the closer we get to the pole, making the density as a function of latitude drop to 0 at the north and south pole (for both the uniform and for the normal case).
The density for $\phi$ will then be the integral of that bivariate density over $\theta$.
Looking at the marginal for $\theta$ (i.e. strips running down at fixed values of $\theta$) makes for four peaks in the density of $\theta$ as you note - indeed this follows directly from the 2D case, but as we now see, it also makes for a pair of peaks in the density of $\phi$ away from the equator, corresponding to a region on the surface of the unit sphere where the corners and upper/lower edges of the cube project.