The basic reason that it is justified is that the time integral is a linear operator.
So let's step way back into the abstract, what you really have is some configuration space $\mathcal C$ of possible configurations of some system, and some space of paths $\mathcal P$ through that, which is some subset of the functions $[0,1]\to \mathcal C$. Inside of this configuration space we may have a time coordinate, position coordinates, velocity or momentum coordinates...
Then we have some functions of these paths, $X:\mathcal P \to \mathbb R.$ Some of them have the format, $$\mathcal A[\mathbf P] = \int_0^1 ds~L\big(\mathbf P(s)\big).$$
We want to do calculus with whole paths, so we want to consider the difference $$\delta \mathcal A = \mathcal A[\mathbf P + \epsilon\mathbf p] - \mathcal A[\mathbf P],~~ \epsilon \approx 0$$ Note that this in theory depends on both the path big-$\mathbf P$ that we use to evaluate the action and the path little-$\mathbf p$ that we use to perturb it, but we're typically going to take the approach of trying to find the big-$\mathbf P$ such that this $\delta \mathcal A$ is zero to first order in $\epsilon$ or so, irrespective of little-$\mathbf p$.
Just putting those two together and using the linearity of the integral, we have $$\delta\mathcal A = \int_0^1 ds~\Big(L\big(\mathbf P(s) + \epsilon \mathbf p(s)\big) - L\big(\mathbf P(s)\big)\Big),$$and the internal term could meaningfully be called $\delta L$. It doesn't mean exactly the same thing as it is applied "pointwise" to the path, but if we were instead to write this as $L[\mathbf P](s)$ or some similar notation then it would indeed be $\delta L[\mathbf P](s)$.
The rest of your expression then follows from the Leibniz property of this $\delta$ operator, which depends on understanding that we are taking a limit etc., which is nontrivial of course.