Consider a liner discrete-time system. Assume we can define it in terms of an input-output relation as follows (you can assume a more general model but it is enough for our purpose):
$$a_0y[n]+a_{1}y[n-1]+\cdots+a_{N}y[n-N]=b_0x[n]+b_{1}x[n-1]+\cdots+b_{M}x[n-M]\tag{1}$$
When the coefficients $\{a_i\}$ and $\{b_i\}$ are constant, we call it a finite-order constant-coefficient ordinary difference equation.
It explains the current output $y[n]$ in terms of a weighted sum of the current and past inputs and the past outputs:
$$y[n]=\frac{-a_{1}y[n-1]-\cdots-a_{N}y[n-N]+b_0x[n]+b_{1}x[n-1]+\cdots+b_{M}x[n-M]}{a_0}$$
It is very similar to a differential equation in continuous time.
Dealing with such equation to find the output of the system can become complicated. It is a recurrence relation. It defines the values recursively and for arbitrary inputs it is not straightforward to express the output in a closed-form representation.
How to deal with it in an easier way?
Consider the following transform (assume the sum exists):
$$\mathcal{Z}(x[n])=\sum_{n=-\infty}^{+\infty}x[n]z^{-n}$$
You are right, it is just a transform that accept $x[n]$ and gives you $X(z)$.
But it has a useful property. Lets calculate the $z$-transform of $x[n-\alpha]$:
$$\begin{align}\mathcal{Z}(x[n-\alpha])&=\sum_{n=-\infty}^{+\infty}x[n-\alpha]z^{-n}\\[10pt]
&=\sum_{n'=-\infty}^{+\infty}x[n']z^{-(n'+\alpha)}\tag{2}\\[10pt]
&=\left(\sum_{n'=-\infty}^{+\infty}x[n']z^{-n'}\right)z^{-\alpha}\\
&=\left(\mathcal{Z}(x[n])\right)z^{-\alpha}
\end{align}$$
where in $(2)$ I changed the variable $n'=n-\alpha \Rightarrow n=n'+\alpha$.
Assume $\mathcal{Z}(x[n])=X(z)$. We have seen that a property of $z$-transform is: $$\boxed{\mathcal{Z}(x[n-\alpha])=z^{-\alpha}X(z)}$$
Applying this to the difference equation $(1)$ makes it (note that it is a liner transform and can be applied term by term):
$$a_0Y(z)+a_1z^{-1}Y(z)+\cdots+a_Nz^{-N}Y(z)=b_0X(z)+b_1z^{-1}X(z)+\cdots+b_Mz^{-M}X(z)$$
Note that all terms are expressed in terms of $Y(z)$ and $X(z)$. So we can say that
$z$-transform reduces usually difficult-to-handle recurrence relations to much easier to manipulate algebraic relations.
With an algebraic relation, we can factor the terms:
$$Y(z)\left(a_0+a_1z^{-1}+\cdots+a_Nz^{-N}\right)=X(z)\left(b_0+b_1z^{-1}+\cdots+b_Mz^{-M}\right)$$
and consequently,
$$Y(z)=X(z)\frac{b_0+b_1z^{-1}+\cdots+b_Mz^{-M}}{a_0+a_1z^{-1}+\cdots+a_Nz^{-N}}$$
where $$H(z)\triangleq\frac{b_0+b_1z^{-1}+\cdots+b_Mz^{-M}}{a_0+a_1z^{-1}+\cdots+a_Nz^{-N}}$$
is called the system function. So we distilled all the system into an algebraic expression as $H(z)$. Hence, the zeros and poles of $H(z)$ have a direct impact on the output.
There are also other reasons that make $z$-transform and $H(z)$ important, and we would prefer working with the $z$-domain rather than time-domain.
Assume the input to the system is $x[n]=z_{k}^n$, where $z_k$ is just a complex number (don't mix it with the $z$ variable in the $z$-transform at this moment! Also there are some restrictions on $z_k$ as we will see next).
The output of the LTI system can be calculated by convolution of the input and the impulse response $h[n]$:
$$\begin{align}
y[n]&=h[n]*x[n]\\
&=\sum_{m=-\infty}^{+\infty} h[m] x[n-m]\\
&=\sum_{m=-\infty}^{+\infty} h[m] z_k^{n-m}\\
&=z_k^n \sum_{m=-\infty}^{+\infty} h[m] z_k^{-m}\\
&=z_k^n H(z_k)
\end{align}$$
What does it really mean?
The output of the system is an scaled version of the input. So the signal $z_k^n$ passes through the system and the only impact on it is a scaling. Note that the scaling factor is the system function when evaluated at $z_k$. Why is it important?
If we could decompose an arbitrary signal (not of the form $z_k^n$) into several components of the form $z_k^{n}$, then we can easily calculate the output of the system to each of the components, and then add them together to find the overall output (benefiting from the superposition property of LTI systems).
How can we decompose the signal $x[n]$ into components of the form $z_k^n$?
It is easy as the following (assume the sum exists for now)
$$\left(\sum_{n=-\infty}^{+\infty}x[n]z_k^{-n}\right)z_k^n$$
Let's refer to this component as $X(z_k)z_k^n$. So the output to this component is $y_k[n]=H(z_k)X(z_k)z_k^n$. We can have as many components as possible. Potentially, we have thousands of $z_k$. So we refer to all of them by a variable $z$ and the component-wise output becomes $H(z)X(z)z^n$. The overall output $y[n]$ is in form of an integral (since we need to add them up). It can be shown that $y[n]$ is given by the following contour integral:
$$y[n]=\frac{1}{2\pi j}\oint_C\frac{H(z)X(z)z^n}{z}dz$$
where $C$ is a counterclockwise path inside the ROC (where $X(z)$ and $H(z)$ exist) that encircles the origin, which is the same as the inverse $z$-transform of $X(z)H(z)=Y(z)$.
Regarding your second question, in the definition of $z$-transform (when it exists on the complex unit circle) if you choose $z=e^{j\omega}$ it becomes the definition of the Discrete-Time Fourier transform Hence is the frequency-domain interpretation.
