Why is the digital frequency response taken on the unit circle, while the analog is taken along the imaginary axis?

One way to see this is to consider the Laplace transform of a sampled signal:

$$x_d(t)=\sum_{n=0}^{\infty}x(nT)\delta(t-nT)\tag{1}$$

where I've assumed that $x(t)$ starts at $t=0$. $T$ is the sampling period, and $\delta(t)$ is the Dirac delta impulse.

Taking the Laplace transform of $(1)$ gives

\begin{align*} X_d(s) &= \int_0^{\infty}\sum_{n=0}^{\infty}x(nT)\delta(t-nT)e^{-st}dt \\ &= \sum_{n=0}^{\infty}x(nT)\int_{0}^{\infty}\delta(t-nT)e^{-st}dt \\ &= \sum_{n=0}^{\infty}x(nT)e^{-snT}\tag{2} \end{align*}

Evaluating $(2)$ on the imaginary axis $s=j\omega$ results in the DTFT of the sampled signal:

$$X_d(j\omega)=\sum_{n=0}^{\infty}x(nT)e^{-j\omega nT}\tag{3}$$

In $(3)$ it is obvious that by varying $\omega$ we move along the unit circle. $X_d(j\omega)$ is periodic with period $2\pi/T$, and for this reason it is frequently written as a function of $e^{j\omega}$.

Setting $e^{sT}=z$ in $(2)$ gives us the $\mathcal{Z}$-transform of the sequence $x(nT)$.