Method One:
\begin{eqnarray*}
& & \int\frac{d^{4}k}{(2\pi)^{4}}\frac{i}{k^{2}+i\epsilon}e^{-ik\cdot x}\\
& = & \frac{i}{(2\pi)^{4}}\int d^{3}ke^{i\mathbf{k}\cdot\mathbf{x}}\int dk_{0}e^{-ik_{0}x_{0}}\frac{1}{[k_{0}+(|\mathbf{k}|-i\epsilon)][k_{0}-(|\mathbf{k}|-i\epsilon)]}\\
& = & \frac{i}{(2\pi)^{4}}\int d^{3}ke^{i\mathbf{k}\cdot\mathbf{x}}\bigg(\theta(x_{0})(-2\pi i)\frac{1}{2|\mathbf{k}|}e^{-i|\mathbf{k}|x_{0}}+\theta(-x_{0})(2\pi i)\frac{1}{-2|\mathbf{k}|}e^{i|\mathbf{k}|x_{0}}\bigg)\\
& = & \frac{\pi}{(2\pi)^{4}}\int d^{3}k\frac{1}{|\mathbf{k}|}e^{i\mathbf{k}\cdot\mathbf{x}}e^{-i|\mathbf{k}||x_{0}|}\\
& = & \frac{\pi}{(2\pi)^{4}}2\pi\int_{0}^{\infty}d|\mathbf{k}||\mathbf{k}|^{2}\frac{1}{|\mathbf{k}|}e^{-i|\mathbf{k}||x_{0}|}\int_{-1}^{1}dye^{i|\mathbf{k}||\mathbf{x}|y}\\
& = & \frac{\pi}{(2\pi)^{4}}2\pi\int_{0}^{\infty}d|\mathbf{k}||\mathbf{k}|e^{-i|\mathbf{k}||x_{0}|}\frac{1}{i|\mathbf{k}||\mathbf{x}|}(e^{i|\mathbf{k}||\mathbf{x}|}-e^{-i|\mathbf{k}||\mathbf{x}|})\\
& = & \frac{\pi}{(2\pi)^{4}}2\pi\frac{1}{i|\mathbf{x}|}\int_{0}^{\infty}d|\mathbf{k}|(e^{-i|\mathbf{k}|(|x_{0}|-|\mathbf{x}|)}-e^{-i|\mathbf{k}|(|x_{0}|+|\mathbf{x}|)})\\
& = & \frac{\pi}{(2\pi)^{4}}2\pi\frac{1}{i|\mathbf{x}|}\bigg[\bigg(\pi\delta(|x_{0}|-|\mathbf{x}|)-\mathscr{P}\frac{i}{|x_{0}|-|\mathbf{x}|}\bigg)-\bigg(\pi\delta(|x_{0}|+|\mathbf{x}|)-\mathscr{P}\frac{i}{|x_{0}|+|\mathbf{x}|}]\bigg)\bigg]\\
& = & \frac{\pi}{(2\pi)^{4}}2\pi\frac{1}{i|\mathbf{x}|}\bigg[\pi\delta(|x_{0}|-|\mathbf{x}|)+\mathscr{P}\bigg(\frac{i}{|x_{0}|+|\mathbf{x}|}-\frac{i}{|x_{0}|-|\mathbf{x}|}\bigg)\bigg]\\
& = & \frac{\pi}{(2\pi)^{4}}2\pi\frac{1}{i|\mathbf{x}|}\bigg(2\pi|\mathbf{x}|\delta(x^{2})-2i|\mathbf{x}|\mathscr{P}\frac{1}{x^{2}}\bigg)\\
& = & \frac{\pi}{(2\pi)^{4}}2\pi\frac{1}{i|\mathbf{x}|}(-2i|\mathbf{x}|)\bigg(\mathscr{P}\frac{1}{x^{2}}+i\pi\delta(x^{2})\bigg)\\
& = & -\frac{1}{4\pi^{2}}\frac{1}{x^{2}-i\epsilon}
\end{eqnarray*}
i.e.
\begin{eqnarray*}
\int\frac{d^{4}k}{(2\pi)^{4}}\frac{i}{k^{2}+i\epsilon}e^{-ik\cdot x} & = & -\frac{1}{(2\pi)^{2}}\frac{1}{x^{2}-i\epsilon}
\end{eqnarray*}
where we have used
\begin{eqnarray*}
\frac{1}{x+i\epsilon} & = & \mathscr{P}\frac{1}{x}-i\pi\delta(x)
\end{eqnarray*}
\begin{eqnarray*}
\int_{0}^{\infty}e^{ikx}dk & = & \lim_{\epsilon\to0^{+}}\int_{0}^{\infty}e^{ik(x+i\epsilon)}dk=\lim_{\epsilon\to0^{+}}\frac{i}{x+i\epsilon}=\pi\delta(x)+\mathscr{P}\frac{i}{x}
\end{eqnarray*}
\begin{eqnarray*}
\delta(x^{2}) & = & \delta(x_{0}^{2}-\mathbf{x}^{2})=\frac{1}{2|\mathbf{x}|}[\delta(x_{0}-|\mathbf{x}|)+\delta(x_{0}+|\mathbf{x}|)]\\
& = & \frac{1}{2|\mathbf{x}|}[\theta(x_{0})\delta(x_{0}-|\mathbf{x}|)+\theta(-x_{0})\delta(x_{0}+|\mathbf{x}|)]\\
& = & \frac{1}{2|\mathbf{x}|}\delta(|x_{0}|-|\mathbf{x}|)
\end{eqnarray*}
Method Two:
We can also see that
\begin{eqnarray*}
& & \int\frac{d^{4}k}{(2\pi)^{4}}\frac{i}{k^{2}+m^{2}+i\epsilon}e^{-ik\cdot x}\\
& = & \frac{i}{(2\pi)^{4}}\int d^{3}ke^{i\mathbf{k}\cdot\mathbf{x}}\int dk_{0}e^{-ik_{0}x_{0}}\frac{1}{[k_{0}+(E_{k}-i\epsilon)][k_{0}-(E_{k}-i\epsilon)]}\\
& = & \frac{i}{(2\pi)^{4}}\int d^{3}ke^{i\mathbf{k}\cdot\mathbf{x}}\bigg(\theta(x_{0})(-2\pi i)\frac{1}{2E_{k}}e^{-iE_{k}x_{0}}+\theta(-x_{0})(2\pi i)\frac{1}{-2E_{k}}e^{iE_{k}x_{0}}\bigg)\\
& = & \frac{i(-2\pi i)}{2(2\pi)^{4}}\int d^{3}k\frac{1}{E_{k}}e^{i\mathbf{k}\cdot\mathbf{x}}e^{-iE_{k}|x_{0}|}\\
& = & \frac{i(-2\pi i)(2\pi)}{2(2\pi)^{4}}\int_{0}^{\infty}dkk^{2}\frac{1}{E_{k}}e^{-iE_{k}|x_{0}|}\int_{-1}^{1}dye^{ik|\mathbf{x}|y}\\
& = & \frac{i(-2\pi i)(2\pi)}{2(2\pi)^{4}}\int_{0}^{\infty}dkk^{2}\frac{1}{E_{k}}e^{-iE_{k}|x_{0}|}\frac{1}{ik|\mathbf{x}|}(e^{ik|\mathbf{x}|}-e^{-ik|\mathbf{x}|})\\
& = & \frac{i(-2\pi i)(2\pi)2}{2(2\pi)^{4}}\int_{0}^{\infty}dkk\frac{1}{E_{k}}e^{-iE_{k}|x_{0}|}\frac{\mathrm{sin}(k|\mathbf{x}|)}{|\mathbf{x}|}\\
& = & \frac{1}{(2\pi)^{2}}\int_{0}^{\infty}dkk\frac{1}{E_{k}}e^{-iE_{k}|x_{0}|}\frac{\mathrm{sin}(k|\mathbf{x}|)}{|\mathbf{x}|}\\
& \equiv & \frac{1}{(2\pi)^{2}}\bigg(\theta(x^{2})\times\mathrm{I}+\theta(-x^{2})\times\mathrm{II}\bigg)
\end{eqnarray*}
If $x^{2}>0$, we can choose a frame with $x^{\mu}=(x_{0},\mathbf{0})$
and $x^{2}=x_{0}^{2}$. We can see
\begin{eqnarray*}
\mathrm{I} & = & \int_{0}^{\infty}dkk\frac{1}{E_{k}}e^{-iE_{k}|x_{0}|}\frac{\mathrm{sin}(k|\mathbf{x}|)}{|\mathbf{x}|}=\int_{0}^{\infty}dkk^{2}\frac{1}{E_{k}}e^{-iE_{k}|x_{0}|}\\
& = & \int_{m}^{\infty}dE_{k}\sqrt{E_{k}^{2}-m^{2}}e^{-iE_{k}|x_{0}|}\\
& = & m^{2}\int_{1}^{\infty}dt\sqrt{t^{2}-1}e^{-im|x_{0}|t},\ \bigg[\text{note: }a\equiv m|x_{0}|=m\sqrt{x^{2}}\bigg]\\
& = & m^{2}\int_{1}^{\infty}dt(t^{2}-1)\frac{e^{-iat}}{\sqrt{t^{2}-1}}\\
& = & -m^{2}(\frac{\partial^{2}}{\partial a^{2}}+1)\int_{1}^{\infty}dt\frac{e^{-iat}}{\sqrt{t^{2}-1}}\\
& = & \frac{i\pi m^{2}}{2}[H_{0}^{(2)\prime\prime}(a)+H_{0}^{(2)}(a)]
\end{eqnarray*}
where we have used
\begin{eqnarray*}
\int_{1}^{\infty}dt\frac{e^{-iat}}{\sqrt{t^{2}-1}} & = & -\frac{i\pi}{2}J_{0}(a)-\frac{\pi}{2}N_{0}(a)=-\frac{i\pi}{2}H_{0}^{(2)}(a),\ (a>0)
\end{eqnarray*}
From
\begin{eqnarray*}
Z_{\nu}^{\prime} & = & Z_{\nu-1}-\frac{\nu}{x}Z_{\nu}\\
Z_{\nu}^{\prime} & = & -Z_{\nu+1}+\frac{\nu}{x}Z_{\nu}
\end{eqnarray*}
we can see
\begin{eqnarray*}
Z_{0}^{\prime} & = & -Z_{1}\\
Z_{0}^{\prime\prime} & = & -Z_{1}^{\prime}=-(Z_{0}-\frac{1}{x}Z_{1})
\end{eqnarray*}
i.e.
\begin{eqnarray*}
Z_{0}^{\prime\prime}+Z_{0} & = & \frac{1}{x}Z_{1}
\end{eqnarray*}
So we have
\begin{eqnarray*}
\mathrm{I} & = & \frac{i\pi m^{2}}{2}[H_{0}^{(2)\prime\prime}(a)+H_{0}^{(2)}(a)]=\frac{i\pi m^{2}}{2a}H_{1}^{(2)}(a)=\frac{i\pi m}{2\sqrt{x^{2}}}H_{1}^{(2)}(m\sqrt{x^{2}})
\end{eqnarray*}
If $x^{2}<0$, we can choose a frame with $x^{\mu}=(0,\mathbf{x})$
and $x^{2}=-|\mathbf{x}|^{2}$. We can see
\begin{eqnarray*}
\mathrm{II} & = & \int_{0}^{\infty}dkk\frac{1}{E_{k}}e^{-iE_{k}|x_{0}|}\frac{\mathrm{sin}(k|\mathbf{x}|)}{|\mathbf{x}|}=\frac{1}{|\mathbf{x}|}\int_{0}^{\infty}dkk\frac{\mathrm{sin}(k|\mathbf{x}|)}{\sqrt{m^{2}+k^{2}}}\\
& = & \frac{m}{|\mathbf{x}|}\int_{0}^{\infty}dt\frac{t\mathrm{sin}(m|\mathbf{x}|t)}{\sqrt{1+t^{2}}},\ \bigg[\text{note: }b\equiv m|\mathbf{x}|=m\sqrt{-x^{2}}\bigg]\\
& = & \frac{m^{2}}{b}\int_{0}^{\infty}dt\frac{t\mathrm{sin}(bt)}{\sqrt{1+t^{2}}}\\
& = & -\frac{m^{2}}{b}\frac{\partial}{\partial b}\int_{0}^{\infty}dt\frac{\mathrm{cos}(bt)}{\sqrt{1+t^{2}}}\\
& = & -\frac{m^{2}}{b}K_{0}^{\prime}(b)
\end{eqnarray*}
where we have used
\begin{eqnarray*}
\int_{0}^{\infty}dt\frac{\mathrm{cos}(bt)}{\sqrt{1+t^{2}}} & = & K_{0}(b),\ (x>0)
\end{eqnarray*}
From
\begin{eqnarray*}
Z_{\nu}^{\prime} & = & Z_{\nu-1}-\frac{\nu}{x}Z_{\nu}\\
Z_{\nu}^{\prime} & = & -Z_{\nu+1}+\frac{\nu}{x}Z_{\nu}
\end{eqnarray*}
we can see
\begin{eqnarray*}
Z_{0}^{\prime} & = & -Z_{1}
\end{eqnarray*}
So we get
\begin{eqnarray*}
\mathrm{II} & = & -\frac{m^{2}}{b}K_{0}^{\prime}(b)=\frac{m^{2}}{b}K_{1}(b)=\frac{m}{\sqrt{-x^{2}}}K_{1}(m\sqrt{-x^{2}})
\end{eqnarray*}
Finally we have
\begin{eqnarray*}
\int\frac{d^{4}k}{(2\pi)^{4}}\frac{i}{k^{2}+m^{2}+i\epsilon}e^{-ik\cdot x} & = & \frac{1}{(2\pi)^{2}}\bigg(\theta(x^{2})\times\frac{i\pi m}{2\sqrt{x^{2}}}H_{1}^{(2)}(m\sqrt{x^{2}})+\theta(-x^{2})\times\frac{m}{\sqrt{-x^{2}}}K_{1}(m\sqrt{-x^{2}})\bigg)\\
& = & \frac{m}{(2\pi)^{2}\sqrt{|x^{2}|}}\bigg(\theta(x^{2})\frac{i\pi}{2}H_{1}^{(2)}(m\sqrt{x^{2}})+\theta(-x^{2})K_{1}(m\sqrt{-x^{2}})\bigg)
\end{eqnarray*}
then we have
\begin{eqnarray*}
& & \int\frac{d^{4}k}{(2\pi)^{4}}\frac{i}{k^{2}+i\epsilon}e^{-ik\cdot x}\\
& = & \lim_{m\to0}\int\frac{d^{4}k}{(2\pi)^{4}}\frac{i}{k^{2}+m^{2}+i\epsilon}e^{-ik\cdot x}\\
& = & \lim_{m\to0}\frac{m}{(2\pi)^{2}\sqrt{|x^{2}|}}\bigg(\theta(x^{2})\frac{i\pi}{2}\frac{2i}{\pi m\sqrt{x^{2}}}+\theta(-x^{2})\frac{1}{m\sqrt{-x^{2}}}\bigg)\\
& = & \frac{1}{(2\pi)^{2}}\bigg(-\theta(x^{2})\frac{1}{x^{2}}+\theta(-x^{2})\frac{1}{-x^{2}}\bigg)\\
& = & -\frac{1}{(2\pi)^{2}}\frac{1}{x^{2}}
\end{eqnarray*}
where we have used
\begin{eqnarray*}
H_{1}^{(2)}(x) & = & J_{1}(x)-iN_{1}(x)\rightarrow\frac{2i}{\pi x}\ \ \text{as}\ \ x\rightarrow0\\
K_{1}(x) & \rightarrow & \frac{1}{x}\ \ \text{as}\ \ x\rightarrow0
\end{eqnarray*}
So we get
\begin{eqnarray*}
\int\frac{d^{4}k}{(2\pi)^{4}}\frac{i}{k^{2}+i\epsilon}e^{-ik\cdot x} & = & -\frac{1}{(2\pi)^{2}}\frac{1}{x^{2}}
\end{eqnarray*}
with $\frac{1}{x^{2}-i\epsilon}$ replaced by $\frac{1}{x^{2}}$.
In short, yes there's an extra minus sign when we reorder the terms. I've just encountered this myself and I found going very slowly helps to understand where the minus sign comes in. Let's use index notation - for some anticommuting $\psi$ (column) and $\bar{\psi}$ (row), and some matrix $M$, consider
$$\begin{align}
\bar{\psi} M \psi &= \bar{\psi}_i M_{ij} \psi_j\\
&= \bar{\psi}_i \psi_j M_{ij} &&\text{rearranging terms}\\
&= -\psi_j \bar{\psi}_i M_{ij} &&\text{anticommuting } \bar{\psi}, \psi \tag{!}\\
&= -\psi_i \bar{\psi}_j M_{ji} &&\text{relabelling } i \leftrightarrow j\\
&= - \psi_i M_{ji} \bar{\psi}_j &&\text{another rearrangment}\\
&= - \psi_i (M^T)_{ij} \bar{\psi}_j &&\text{definition of transpose}\\
&= -\psi^T M^T \bar{\psi}^T &&`\text{row $\times$ matrix $\times$ column'}
\end{align}$$
which reproduces the usual 'transposes in inverse order' but with an extra minus sign because of the anticommuting nature of the $\bar \psi$ and $\psi$ at $(!)$.
When it comes to something being equal to its own transpose, that's still correct - however, we don't have the usual 'transpose of a product' formula and instead we get extra minus signs, e.g.
$$\bar{\psi} \psi = (\bar \psi \psi)^T = -\psi^T \bar{\psi}^T$$
Best Answer
The free scalar and fermion propagator is $$ G_\psi(x,y) = \int \frac{d^dp}{(2\pi)^d} \frac{-i(\gamma^\mu p_\mu + m)}{ p^2 + m^2 - i \epsilon} e^{- i p \cdot ( x - y ) } $$ The scalar propagator is $$ G_\phi(x,y) = \int \frac{d^dp}{(2\pi)^d} \frac{-i}{ p^2 + m^2 - i \epsilon} e^{- i p \cdot ( x - y ) } $$ Clearly, $$ G_\psi(x,y) = ( i \gamma^\mu \partial_\mu -m)G_\phi(x,y)~. $$
PS - In any dimension, the gamma matrices are defined to satisfy $\{ \gamma^\mu , \gamma^\nu \} = - 2 \eta^{\mu\nu}$.
PPS - I am using metric signature $(-,+,+,+,\cdots)$ in this answer.