How/why are the $\mathcal Z$-transform and unit delays related?
Mathematically we can easily show that the z-transform of a unit cycle delay is $z^{-1}$, Just like the s-transform of a time delay $\tau$ is $e^{-\tau s}$. I would like to add some additional insights that may lead to a more intuitive understanding:
Why z??
NOTE: To really understand how the z-transform and unit delays are related (intuitively), it really helps to understand the Laplace domain, what it is used for, and how it maps to the z-domain. For brevity I am including the first graphic covering a main point needed and the following paragraph giving a high level perspective, without going into further details of how we use either plane, but the more you know about that, the more the rest will make sense of what I will show after that.-- For further details on s and z along with a great graphic by Julius Smith that inspired my plots below, see  https://ccrma.stanford.edu/~jos/mdft/Comparing_Analog_Digital_Complex.html)
To help understand $z^{-1}$ it is useful to understand why we use the transform into the s and z plane. This is an entire topic on its own but I will at least make note of this very interesting observation: The process of multiplying (complex conjugate multiplication to be specific) and then integrating the product for two functions in the continuous time domain, or multiplying and then adding the product in the discrete time domain can be viewed as &quot;correlation&quot;. (In actuality it is impressive how many things we know in DSP are a form of this description of correlation!).  Take a close look at the equations for the Fourier Series Expansion (FSE), the Fourier Transform (FT), Laplace Transform (LT) and Z-Tranform (ZT) equations and consider how we are doing correlations to each possible frequency for the FSE ($e^{jn\omega_nt}$), each point on the frequency axis for the FT ($e^{j\omega t}$), and each point on the complex surface for S ($e^{s t}$) and Z ($z^n$):
Fourier Series Coefficients
(scaled correlation of x(t) with $e^{jn\omega_0t}$, x(t) is only defined from 0 to T):
$$c_n = \frac{1}{T}\int_0^T x(t)e^{-jn\omega_0t}dt $$
Fourier Transform (correlation of x(t) with $e^{j\omega t}$) :
$$X(\omega) = \int_{t=0}^\infty x(t)e^{-j\omega t}dt$$
Laplace Transform (correlation of x(t) with $e^{st}$):
$$X(s) = \int_{t=0}^\infty x(t)e^{-st}dt$$
Z-Transform (correlation of x[n] with $z^n$):
$$X(z) = \sum_{n=0}^\infty x[n]z^{-n}$$
(x(t) and x[n] assumed to be causal, hence the integrations and summations start at 0).
(Note that for the Laplace and Z Transform it is not strictly a correlation the way I defined it since s and z have real components, so $e^{st}$ is not simply the complex conjugate of $e^{-st}$; both the real and imaginary exponent signs change.)
In each case, we have mapped our time domain function by describing it in terms of other variables, with the primary benefit of revealing its extrema (the poles and zeros) that can uniquely and compactly describe a (linear time invariant) system, opening up a whole world of elegant problem solving (notably solving integro-differential equations with simple algebra, and describing/ manipulating systems based on their pole and zero locations).
So we convert our discrete time domain sequence into the continuous z-domain by correlating it with every possible value of $z^n$ where z is any complex number, thereby decomposing our waveform into an infinite number of z's each with its own complex weight (magnitude and phase). This is the surface of our function on the z domain, which we usually just view the pole and zero locations (as that is sufficient to determine the magnitude and phase of every other point on that surface!).
The unit circle in the z-plane is the frequency axis:
Below demonstrates one way, out of many, to map from the s-plane to the z-plane. For purposes here, the only point I want to make is the reminder that the vertical axis in s which represents frequency $\omega$ (the Fourier transform is the function H(s) when you set $s=j\omega$),  maps to the unit circle $|z|=1$ on the z-plane.
(This mapping is specifically the matched-z transform, other common approaches are the bilinear transform, which I use most often, and the method of impulse invariance.  In all these mappings, what is consistent is the unit circle in z is the frequency axis (imaginary axis in s) and the left half plane in s maps to the inside of the unit circle in z. The answer on this link was very good for more details into the mapping techniques and why: Are there alternatives to the bilinear transform?)
Additionally this answer by Matt L describes additional mapping techniques where the $j\omega$ axis in the s-plane does NOT map to the unit circle in z-plane: z-Transform Methods: Definition vs. Integration Rule




Time Delay of One Sample = $z^{-1}$:
The math to show that the z-transform of a unit cylce delay is $z^{-1}$ is extremely simple (which is why we love the z-transform; we could certainly compute everything in discrete time systems using Laplace if we wanted the punishment).
$X(z) = Z\left\{x(n)\right\} =\sum_{n=0}^\infty x(n)z^{-n} $
$Z\left\{x(n-1)\right\} = X(z)z^{-1}$
So knowing this I provide the following graphics below to make the main point: as we increase powers of z, we are moving forward in time. The first graphic is of a shift register, where we often see each register denoted as $z^{-1}$, indicating that that element provides a 1 sample delay. Therefore the contents of the shift register at any given moment contain the current sample in the left most position, and all the previous samples as we move forward to the right (so we are moving backward in time reviewing the contents of our waveform in the shift register at that moment).  The second graphic shows the same waveform, but as we would be used to displaying it, on an axis that is moving forward in time, and that axis has increasing powers of z.





What is z???
Next, let's talk about z and what it is. For purposes here,the important point to know is that z is any complex number.



Now, knowing that z can be any complex number, and increasing powers of z represents samples as we move forward in time (and with the reminder when you multiply two complex numbers, you add the phases and multiply the magnitudes), let's look at $z^n$ as $n=0,1,2,3...$ for various values of z :








Notice in these plots how when z as any complex number is limited to the unit circle, the waveform that results for $z{^n}$ is a constant amplitude complex frequency, just as we introduced in the first plot how the unit circle represents the frequency axis.  When z=1, the frequency is DC, and when z is magnitude 1 with a very small angle, the frequency is small, and as the angle of z increases, the frequency increases. A delay in time of a complex number with magnitude 1 results in a rotation (the angles add), so if we continuously add delay with time, we will get a continuous rotation (frequency is change in phase vs change in time). This hopefully will give more insight in tying together an interpretive meaning for z, and $z^{-1}$ and how it relates to a unit sample delay.