Main idea
Indeed as @JeanbaptisteRoux pointed out in a comment, the motivation is that you can think of the determinant of an operator to a positive power as coming out of a fermionic integral. To illustrate this without the necessity of all the technology of non-abelian gauge theories let me give you a toy example that could be treated with Faddeev-Popov ghosts; a massless particle in one-dimension.
Toy example
Consider then a free massless field:
$$\newcommand{\d}{\mathrm{d}}\newcommand{\pd}{\partial}\newcommand{\D}{\mathrm{D}}S[\phi] := \frac12 \int_0^1\d{\tau}\ \big(\pd_\tau\phi(\tau)\big)^2.$$
For concreteness in the calculations I will be in the Euclidean signature, so this is the Euclidean action and $\tau$ is Euclidean time.
This action has a shift symmetry $\phi\mapsto \phi+\lambda$ where $\alpha$ is a constant. You can gauge the symmetry by considering a gauge field $A_\tau$ that transforms as $A_\tau\mapsto A_\tau^\lambda :=A_\tau + \pd_\tau \lambda$ under the symmetry, so the gauged action is
$$S[\phi, A] = \frac12\int_0^1 \d{\tau}\ \big(\pd_\tau \phi(\tau) - A_\tau\big)^2.$$
To compute the partition function of the theory you need to integrate over $\phi$ and $A$ and divide by the volume of the gauge group — in this case this is $\mathcal{G}=\operatorname{Maps}(\mathbb{S}^1\to\mathbb{R})$ — to avoid overcounting gauge configurations. All in all you have
$$ Z_\text{gauge} = \int \frac{\D{\phi}\D{A}}{\operatorname{vol}(\mathcal{G})}\ \exp(-S[\phi,A]).$$
Now comes the first Faddeev-Popov trick. Choose a specific gauge-field that you like, $\hat{A}_\tau$ (this could be for example $\hat{A}_\tau=0$ in this case) and note that
$$\DeclareMathOperator{\det}{det} \det(\pd_\tau) \int\D{\lambda}\ \delta\!\left(A-\hat{A}^\lambda\right) = 1, \tag{FP1}$$
where $\delta(\cdot)$ is a functional generalisation of the delta-function. To convince yourself that this is true, remember the property of the Dirac delta-function
$$\delta\big(f(x)) = \left\vert\frac{1}{f'(x_0)}\right\vert\delta(x-x_0),$$
where $x_0$ is the only zero of $f(x)$; then generalise for a multidimensional delta-function, to finally generalise for a functional delta-function and land on (FP1). In the above I have been very cavalier with zero-modes; the correct equation needs to deal with zero-modes separately, but that is a story for another day (question).
Cool, now insert $1$ in the guise of (FP1) in your partition function to get
$$Z_\text{gauge} \overset{\text{(FP1)}}{=\!\!=\!\!=} \int \frac{\D{\phi}\D{A}}{\operatorname{vol}(\mathcal{G})}\ \det(\pd_\tau) \int\D{\lambda}\ \delta\!\left(A-\hat{A}^\lambda\right)\exp(-S[\phi,A]).$$
Next thing to notice is that you can do the integral over $A$ (let's choose for simplicity $\hat{A}=0$) because of the delta-function and then everything that remains is gauge-invariant so you can do the $\lambda$-integral too. What you get is
$$\int\D\lambda = \operatorname{vol}(\mathcal{G}),$$
cancelling exactly the volume in the denominator. This was probably a bit too fast but if you take a bit of time at every non-trivial step you can convince yourself it is true.
All in all, what you're left with is
$$Z_\text{gauge} = \det(\pd_\tau)\int\D{\phi} \exp(-S[\phi]).$$
And here comes the second Faddeev-Popov trick. Note the following
$$\det(\pd_\tau) = \int \D{\bar{c}}\D{c}\ \exp\!\left(-\int_0^1\d{\tau}\ \bar{c}(\tau)\pd_\tau c(\tau)\right),\tag{FP2}$$
where $c$ is a complex Grassmann-valued field; a ghost; a spin-0 fermion. This is simply a property of Grasmannian integration. Again, I have sacrificed (important) subtleties of zero-modes in order to give the big picture, without many digressions.
So in total you can think of the partition function of your gauge theory as coming from a theory containing a weird bunch of spin-0 fermions, i.e.
\begin{align}Z_\text{gauge} &= \int \D{\phi}\D{\bar{c}}\D{c}\ \exp(-S_\text{tot}[\phi,c]),\qquad \text{where}\\
S_\text{tot}[\phi,c] &= \int_0^1\d{\tau}\ \left[\frac12\big(\pd_\tau\phi(t)\big)^2+\bar{c}(\tau)\pd_\tau c(\tau)\right].
\end{align}
In QCD
The exact same steps lead to the introduction of ghosts in QCD. What differs is that the gauge transformation is, instead that the gluons (the analogues of $A_\tau$) have $\mathrm{SU}(N_\text{colour})$ non-abelian gauge transformations, and that the gauge-fixing condition (the analogue of the choice $\hat{A}=0$) is different. This leads to more complicated descriptions for the total action $S_\text{tot}$, but conceptually it is the same procedure.
Some other comments
To answer some of your other questions:
As you saw above, in the toy example, the description in terms of ghosts is just a choice. It is a convenient choice, because it is the only way to express $\det[\text{operator}]$ in a local way and so it helps with doing perturbative calculations (which is more or less the only thing we can do with interacting theories). We can refrain from using anticommuting variables by simply keeping the determinant as it is, without expressing it locally. The price to pay is that to calculate the determinant, one must solve a, typically very difficult PDE (which is usually impossible).
Theorists aren't afraid of ghosts. They understand where they come from and what their role in life is, and they know the ghosts must die when calculating physical observables.
