This is not a full answer, but I believe that your colleague's result does all the heavy lifting whenever the flow does not have any fixed points. The case near a fixed point of a flow seems more difficult. I'm not sure how that works exactly, but if one can understand it in the $1$-dimensional case (i.e. when the flow has a limit point), this might give some good insight. I broke the answer into several sections to help make it more readable.

$\textbf{Finding the right coordinates}$

In this case, the smooth flow locally foliates your manifold, so around any point $p$ one can use the inverse function theorem to choose local coordinates $\{ x_0, x_1, \ldots, x_{n-1} \}$ for which the flow acts as $\phi_t(x_0) = x_0+t$ and $\phi_t(x_i) = x_i$ otherwise. By Professor Burnhol's argument (which I'm black-boxing for this argument) $f(  x_0, x_1, \ldots, x_{n-1})$ depends smoothly on $x_0$ when all the other coordinates are fixed.

What remains to show is that it depends smoothly on $x_1, \ldots, x_n$. We start with the case where $x_0$ is some fixed value, say 0. In this case, the $A_T f(x)$ is exactly the following:

 $$ \frac{1}{T}\int_0^T f(t, x_1, \ldots, x_n) ~dt$$

For the rest of this, I will use the parameter $t$ for $x_0$ and refer to the rest of the coordinates collectively as $\mathbf{x}$.

We know that $A_T f(x)$ depends smoothly on $\mathbf{x}$ and further that $f$ depends smoothly on $t$. What remains to show is that $f$ depends smoothly on $\mathbf{x}$. We can't take derivatives with respect to $\mathbf{x}$, so we need to estimate some difference quotients. I will be using the Euclidean norm on the $\mathbf{x}$ coordinates for my estimates.

$\textbf{The uniform Lipschitz estimate}$

Take two points $\mathbf{x}^1$ and $\mathbf{x}^2$. Since $A_T f(x)$ is smooth, it is locally Lipschitz, and so we have the following estimate:

$$ \frac{1}{T \| \mathbf{x}^1- \mathbf{x}^2 \| } | \int_0^T f(t, \mathbf{x}^1) - f(t, \mathbf{x}^2)  ~dt | < \mathcal{C} $$

For fixed $\mathbf{x}^1$ and $\mathbf{x}^2$, this estimate is uniform for small $T$ by the smoothness of $f$ in the $t$ variable. This uniformity is crucial. We now pick $T>0$ small enough so that that for $\mathbf{x}^1$ and $\mathbf{x}^2$, we have the following:

 $$ \| f(t, \mathbf{x}^i) -f(s, \mathbf{x}^i)  \|_{C^2} < \epsilon$$
for $0<t,s < T$. This estimate may depend on $\mathbf{x}^1$ and $\mathbf{x}^2$. In particular, we don't expect to be able to do this uniformly in $\mathbf{x}$, but that's okay because we are only doing it for two points.

Then, we have the following estimate.

\begin{eqnarray}
| \frac{1}{T} \int_0^T \frac{ f(t, \mathbf{x}^1) - f(t, \mathbf{x}^2)}{\| \mathbf{x}^1- \mathbf{x}^2 \| }  ~dt | & \geq & | \frac{1}{T} 
\frac{ T(f(0,\mathbf{x}^1) -  f(0,\mathbf{x}^2))}{\| \mathbf{x}^1- \mathbf{x}^2 \|} | \\
& & - \frac{1}{T} \left( \frac{\partial f}{\partial t} (0, \mathbf{x}^1)+ \frac{\partial f}{\partial t} (0, \mathbf{x}^2) + 2\epsilon \right) \frac{T^2}{2} \\
& \geq & | \frac{f(0,\mathbf{x}^1) -  f(0,\mathbf{x}^2)}{\| \mathbf{x}^1- \mathbf{x}^2 \|}| - \left( \frac{\partial f}{\partial t} (0, \mathbf{x}^1)+ \frac{\partial f}{\partial t} (0, \mathbf{x}^2) + 2\epsilon \right) T \\
\end{eqnarray}

Letting $T$ be go to zero using the uniformity of the earlier estimate while still fixing $\mathbf{x}^1$ and $\mathbf{x}^2$ for now, our work shows the following.

$$ | \frac{f(0,\mathbf{x}^1) -  f(0,\mathbf{x}^2)}{\| \mathbf{x}^1- \mathbf{x}^2 \|}| \leq \mathcal{C} $$

Notice that this $\mathcal{C}$ is the same constant as before in the estimate of $A_T f(x)$. Repeating this for arbitrary pairs of $\mathbf{x}$, this gives a uniform Lipschitz estimate on $f$ (at $t=0)$ in terms of the $\mathbf{x}$ coordinates. 

$\textbf{A sketch of how to use induction for higher order regularity}$

In fact, we can repeat essentially the same argument but with three points $\mathbf{x}^1$, $\mathbf{x}^2$ and $ \lambda \mathbf{x}^1 + (1- \lambda) \mathbf{x}^2$ where instead of using the difference quotient, we use the second order difference quotient. This should yield a uniform $C^{1,1}$ estimate on $f$ in terms of the smoothness of $A_T f(x)$. With a uniform $C^{1,1}$ estimate, this shows that the function $f$ was actually differentiable in terms of $\mathbf{x}$ all along. You can induct on all orders to get smoothness of $f$ in $\mathbf{x}$. This shows that $f$ is smooth in both $t$ and $\mathbf{x}$.