This is not a full answer, but I believe that your colleague's result does all the heavy lifting away from the fixed points of the flow. 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.

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.

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.

$\textbf{Update with some corrections}$:

Unfortunately, my original answer had a serious gap. Sorry for the mix-up. The argument I gave only works when we have estimates on $A_T f(x)$ that are uniform in $T$. However, all the hard work was not useless as it made clear the properties that a counterexample would have. Without further ado, here is an example for which $A_T f(x)$ is smooth while $f(x)$ is not.

Let $M= \mathbb{R^2} =(t,x)$ with the flow $\phi_s (t, x) := (t+s, x)$.

Let $\Psi$ be the following function

$$\Psi(x) = \begin{cases} 0 & x\leq 0 \\ e^{\frac{-1}{x^2}} & x > 0 \\ \end{cases}$$

In particular, we want $\Psi$ to be a smooth function which vanishes to all orders at $0$. Let $f$ be the following function:

$$f(t,x)= \begin{cases} e^{\frac{-t^2}{(\Psi(x))^2}} & x >0 \\ 0 & x \leq 0 \end{cases} $$

The function $f$ is not smooth. In fact, it is not even continuous at $(0,0)$. To show that $A_T f(x)$ is smooth for any $T$, we must only show that it depends smoothly on $x$ near $(0,0)$. For fixed $T$>0, we have the following:

$$A_T(f(t_1,x))= \frac{1}{T} \int_{t_1}^{t_1+T} f(t,x) dt < \frac{\sqrt{ \pi}}{T} \Psi(x) = o(x^k) ~ \forall k >0 $$

From this, it follows that $\frac{\partial^k A_T f}{\partial x^k}(t,0) \equiv 0$, so the derivative of $A_T f$ exists to all orders at $(0,0)$ for $T>0$. The smoothness of $A_T f(x)$ elsewhere is immediate because $f$ is smooth away from the origin.

For an example where $f(x)$ is continuous but not smooth, we can consider the function $\hat f(t,x) = x \cdot f(t,x)$. Here, the same analysis holds except that $f$ is continuous but not differentiable at $(0,0)$.

$\textbf{Original answer with the some corrections}$

This is not a full answer, but I believe that your colleague's result does all the heavy lifting away from the fixed points of the flow, at least once you have good estimates on $A_T f(x)$. 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.

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 Burnol'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.

In order for the rest of the argument to work, we must assume that this estimate is uniform in the choice of $T$, $\mathbf{x}^1$ and $\mathbf{x}^2$. This allows us to pick $T>0$ small enough so that that for our fixed $\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$.

The choice of $T$ 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.

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}$.

This is not a full answer, but I believe that your colleague's result does all the heavy lifting away from the fixed points of the flow. 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}$.