Taylor's theorem tells us that to first order
$$ f(x+ \delta x) = f(x) + \delta x f'(x) + O(\delta x^2) $$
Can I simply upgrade $\delta x $ to a function and say
$$ f(x+ g(x)) = f(x) + g(x) f'(x) + O(g^2)?$$
I feel like it should be true and I have seen it being used in equation (2.77) of this. Simply replacing $\delta x$ with a function doesn't seem like a solid proof to me however so I am quite uncomfortable doing this.
I have tried to show this myself using Taylor's theorem for a composite function about a point $x_0$:
$$ f(h(x))=f(h(x_0)) + (x-x_0) f'(h(x_0)) h'(x_0) + \ldots $$
where I take $h(x) = x + g(x)$, this yields
$$ f(x+g(x)) = f(x_0 + g(x_0)) +(x-x_0)f'(x_0+g(x_0))(1 +g'(x_0))+\ldots$$
but I am unsure how to proceed.
In that article, particularly in equation ($2.77$) where they have $\xi^{\mu}(x)$, they're thinking of it as "an infinitesimal change in the coordinates $x^{\mu}$". It might be clearer to think of this as being $\epsilon \cdot \zeta^{\mu}(x)$, for a small real parameter $\epsilon$, and $\zeta^{\mu}$ being any function, and then expand in powers of $\epsilon$ (that's really what was intended).
Now, Taylor's theorem in one-dimension says that \begin{align} f(x+h) &= f(x) + f'(x) \cdot h + o(|h|) \quad \text{as $h \to 0$} \end{align} What this means more explicitly is that there is a "remainder term", which is a function $\rho(\cdot)$ such that \begin{align} f(x+h) &= f(x) + f'(x) \cdot h + \rho(h) \cdot |h| \end{align} and $\lim \limits_{h \to 0}\rho(h) = 0$. Now, I'm guessing you want to consider how to expand something like $f(x + \xi(x))$. As mentioned above, you should think of this as $f(x + \epsilon \cdot \zeta(x))$, for a small real number $\epsilon$, and expand in powers of $\epsilon$ (that's really the definition of the variation). Well, we have \begin{align} f(x + \epsilon \cdot \zeta(x)) &= f(x) + f'(x) \cdot \left( \epsilon \cdot \zeta(x)\right) + \rho(\epsilon \cdot \zeta(x)) \cdot |\epsilon \cdot \zeta(x)| \\ &= f(x) + \epsilon \cdot \underbrace{\zeta(x) f'(x)}_{\text{first order in $\epsilon$}} + \underbrace{\bigg(\epsilon |\zeta(x)| \rho(\epsilon\cdot \zeta(x)) \bigg)}_{= o(\epsilon)} \tag{$*$} \end{align} The last term is $o(\epsilon)$ because \begin{align} \lim_{\epsilon \to 0} \dfrac{\epsilon |\zeta(x)| \rho(\epsilon\cdot \zeta(x))}{\epsilon} = \lim_{\epsilon \to 0}|\zeta(x)| \rho(\epsilon\cdot \zeta(x)) = 0 \end{align} where the last equal sign is because $\rho(h) \to 0$ as $h \to 0$.
Another way to get $(*)$ is to directly apply differentiate $f(x + \epsilon \cdot \zeta(x))$ when regarded as a function of $\epsilon$. In other words, define $F(\epsilon) = f(x + \epsilon \cdot \zeta(x))$. Then, by Taylor's theorem (or in this first order case, simply by definition of differentiability), we have \begin{align} F(\epsilon) &= F(0) + \epsilon \cdot F'(0) + o(|\epsilon|) \quad \text{as $\epsilon \to 0$} \end{align} To calculate $F'(0)$, just use the chain rule; it will be $\zeta(x) \cdot f'(x)$. Hence, \begin{align} f(x + \epsilon \cdot \zeta(x)) &= f(x) + \epsilon \cdot \zeta(x) f'(x) + o(\epsilon), \end{align} exactly as we found in $(*)$.
By the way, to see a similar idea of regarding the smallness as controlled by a parameter $\epsilon$, you may want to take a look at this previous answer of mine.