Why propagator pole is associated to the mass? We say that the pole of the all-orders photon propagator, $$\frac{1}{q^2[1+\Pi(q^2)]}$$
doesn't shift if $\Pi(q^2=0)$ is regular. This amounts to say that the photon remains massless to all orders in perturbation theory. Conversely, the fermion propagator,
$$\frac{i}{\not p -m+\Sigma(p)+i\varepsilon}$$
has $\Sigma(p)$ which shifts the pole hence the mass is corrected and must be renormalized.
Why the propagator pole is associated to the mass?
 A: I believe the simplest way to understand this is via the Kallen-Lehmann spectral representation, which basically teaches you about the spectrum of the theory-in a nonperturbative way-by looking at the analytic structure of the propagator in momentum space.
You begin by considering an interacting theory, with an interacting vacuum $|\Omega\rangle$(related to the free vacuum by the usual trick of taking time evolution to infinity). We have an interacting hamiltonian $H$, and we need a complete set of eigenvectors for this to describe our theory. Note that these will also be simultaneous eigenvectors of the momentum, since $[P,H]=0$(this is a statement about translation invariance of the theory). You can look at Peskin's chapter 7 for a detailed derivation, but the essential result is that the free theory completeness relation gets modified to
$$(1)=|\Omega\rangle\langle\Omega|+\sum_{\lambda_0}\int \frac{d^3p}{(2\pi)^3}\frac{1}{2E_p(\lambda)}|\lambda_p\rangle \langle\lambda_p|$$
Where $|\lambda_p\rangle$ is the state obtained by boosting some zero momentum state $|\lambda_0\rangle$. There can be a lot of such zero momentum states, each will give rise to a set of its own boosted states, and so we sum over $\lambda_0$ in the completeness relation. $E_p(\lambda)=\sqrt{|p|^2+m_\lambda^2}$ is the dispersion relation for $|\lambda_p\rangle$, and it is now clear what the interpretation of $m_\lambda$ is- it is the energy of a ZERO MOMENTUM state $\lambda_0$- this is just we call the mass.
We can now massage this equation-look at Peskin for details, to find the 2 point function-
$$\langle\Omega|\phi(x)\phi(y)|\Omega\rangle=\sum_{\lambda}\int\frac{d^4p}{(2\pi)^4}\frac{1}{p^2-m_\lambda^2} e^{-ip(x-y)}Z$$
Where $Z^2=|\langle\Omega|\phi(0)|\lambda_0\rangle|^2$ is what we interpret as the wavefunction renormalization. It is now clear what the poles of this expression are- they are precisely the masses $m_\lambda$. We can write this in momentum space and the fourier transform-the spectral density-takes the form-
$$\rho(M^2)=Z\rho(M^2-m_{1p}^2)+...$$
where the first pole obviously comes at the 1-partucle states' masses-this is the smallest mass so it comes first. There will be further contributions from bound states, and a cut starting from multiparticle states. The punchline is, the non analyticities are associated with masses in the theory. These non analyticities are just the poles of the propagator.
A: Here is a heuristic argument.

*

*It's a fact that a full connected propagator/2-point correlation function $\tilde{G}_c$ of a Lorentz-invariant theory tend to have a simple pole $\tilde{G}_c \propto \frac{1}{p^2-m^2}$, where for the sake of the argument $m$ is some constant.


*Given that a correlation function/scattering amplitude is a measure of how likely a process occur, it is natural to associate an infinity/a pole with the production of a particle with physical mass $m$.


*See also the Källén–Lehmann spectral representation.


*For unstable/quasi-particles, see also the relativistic Breit–Wigner distribution.
