The integral $$G(x,y)=\int \dfrac{d^4p}{(2\pi)^4}\dfrac{e^{-ip(x-y)}}{p^2+m^2}\tag{1}$$
is undefined because of the fact that the denominator $p^2+m^2$ vanishes when $p$ goes on shell. This corresponds to a whole three-dimensional hyperboloid in the integration domain. This is a special feature of the Lorentzian signature of Minkowski spacetime. Indeed, in Euclidean signature, $p^2+m^2$ can never be zero because $p^2\geq 0$ and then (1) uniquely defines the inverse of the scalar Laplace operator. In Lorentzian signature, we must pick an inverse.
To pick an inverse we introduce a prescription that makes (1) well-defined. Each of the prescriptions defines a different function, with different properties.
First of all, let us write (1) a bit differently by using the definition of $p^2=-(p^0)^2+|\vec{p}|^2$. In that case
$$G(x,y)=\int \dfrac{dp^0d^3\vec{p}}{(2\pi)^4}\dfrac{e^{-ip(x-y)}}{-(p^0)^2+|\vec{p}|^2+m^2}\tag{2}.$$
Now let us imagine that we take the integral over $p^0$ first with $\vec{p}$ held fixed. We can see that the singularity is hit when $(p^0)^2=m^2+|\vec{p}|^2$. This happens at $$p^0 = \pm \sqrt{m^2+|\vec{p}|^2}\tag{3}$$
Since the integral over $p^0$ is over the whole real line, and the two solutions in (3) are both sitting on the line, the integration goes through them as expected. You can circumvent this by defining $G(x,y)$ as
$$\int \dfrac{dp^0d^3\vec{p}}{(2\pi)^4}\dfrac{e^{-ip(x-y)}}{-(p^0\pm i\epsilon)^2+|\vec{p}|^2+m^2},\quad \text{or},\quad \int \dfrac{dp^0d^3\vec{p}}{(2\pi)^4}\dfrac{e^{-ip(x-y)}}{-(p^0)^2+|\vec{p}|^2+m^2\pm i\epsilon} \tag{4}.$$
In the first option you get the singularities at $$p^0=\pm i \epsilon + \sqrt{m^2+|\vec{p}|^2},\quad p^0 = \pm i\epsilon -\sqrt{m^2+|\vec{p}|^2}.\tag{5}$$
The singularities are thus either pushed upward or downwards in the imaginary direction. So suppose we take the sign $+$. The singularities are then pushed upwards. You may close the contour up or down. If you close it upwards you enclose the singularities and get a contribution from the associated residues. If you close it downwards you don't get it. In one case you get a non-zero result, in the other you get zero. So what determines the contour choice?
Well, it is essentially the exponential. Observe you have $e^{-ip(x-y)}$ on the integrand. This is $e^{ip^0(x^0-y^0)}e^{-i\vec{p}\cdot (\vec{x}-\vec{y})}$. Now since we are doing contour integration, we will have both real and imaginary parts of $p^0$ along the contour, so these exponentials are $$e^{i{\rm Re}(p^0)(x^0-y^0)-{\rm Im}(p^0)(x^0-y^0)}e^{-i\vec{p}\cdot (\vec{x}-\vec{y})}\tag{6}.$$
This is the key to the choice of contour. If you close the contour upwards, along the big semicircle you will have ${\rm Im}(p^0)\to \infty$. So if $(x^0-y^0)<0$ you will have one exponential $e^{-{\rm Im}(p^0)(x^0-y^0)}$ which badly diverges, while if $(x^0-y^0)>0$ you have one exponential $e^{-{\rm Im}(p^0)(x^0-y^0)}$ which goes to zero and eliminates the big semicircle contribution. So you must close the contour upwards when $x^0-y^0>0$. By the same argument, you must close the contour downwards when $x^0-y^0<0$.
Now, recall that when you pick the plus sign in (5), the poles are shifted upwards. So you only get a non-zero contribution when the contour is closed upwards, which is the case for $x^0-y^0>0$. So for the plus sign in (5) you get something non-zero for $x^0-y^0>0$ and zero for $x^0-y^0<0$. This is the retarded propagator.
When you pick the minus sign in (5) the opposite happens. The poles are shifted downwards and you only get a non-zero contribution when the contour is closed downwards, which is the case for $x^0-y^0<0$. So for the minus sign, you get a non-vanishing result for $x^0-y^0<0$ and zero for $x^0-y^0>0$. This is the advanced propagator.
Finally, consider the second option in (4). Now the singularities are at $$p^0
= \sqrt{m^2+|\vec{p}|^2\pm i\epsilon},\quad p^0 = - \sqrt{m^2+|\vec{p}|^2\pm i\epsilon}\tag{7}.$$
Since $\epsilon$ will be taken to zero, you can expand the square roots and the singularities lie at $$p^0
= \sqrt{m^2+|\vec{p}|^2}\pm \frac{i\epsilon}{2\sqrt{m^2+|\vec{p}|^2}},\quad p^0 = - \sqrt{m^2+|\vec{p}|^2}\mp \frac{i\epsilon}{2\sqrt{m^2+|\vec{p}|^2}}\tag{8}.$$
Regardless of the choice of sign in the $i\epsilon$ prescription, the important point is that now one pole will be shifted upwards and the other downwards. The contour analysis now shows that in both cases $x^0-y^0>0$ where the contour must be closed upwards and $x^0-y^0<0$ where the contour must be closed downwards you will get contributions. The choice of the sign will determine the time-ordering of the resulting correlator.
So, in summary, (1) by itself is ill-defined because the wave operator in Lorentzian signature does not have a unique inverse. To make (1) well-defined you can choose several possible $i\epsilon$ prescriptions as in (4), which are essentially boundary conditions on causal behavior of the inverse. The choice of prescription then picks an inverse of the wave operator. For each such choice, you may then evaluate the integral over $p^0$ using contours and you discover that that choice of prescription is tied to some causal behavior of $G(x,y)$ because of its dependence on $x^0-y^0$. This is why causality leads to different contours when calculating propagators.