Meirovitch says in his "Principles and Techniques of Vibrations" (1997) on p.85:
In the case of holonomic systems, the variation and integration processes are interchangeable (...)
which means that
$\delta \int_{t_1}^{t_2} L(q, \dot{q}) dt = \int_{t_1}^{t_2} \delta L(q,\dot{q}) dt$
subject to
$f(r_1(q), r_2(q),...,r_N(q),t) = 0.$
where $L$ is the Lagrangian, $q$ are generalized coordinates, $r_i$ are the coordinates of points of mass in the configuration space and $f$ is the constraint.
Can anybody tell me why?