Hopefully this is not a repeated question.
The question comes from a physics background: in quantum field theory, a method to compute quantities of interest is to define the so-called quantum effective action, which measures the response of the quantum system to external sources. It is defined through the path integral, $$iW_{\it eff} = \log\int\mathcal{D}\phi e^{iS[\phi] + \int dx J \phi},$$ where $J$ are the external sources. Then the functional derivatives of $W_{\it eff}$ give expectation values, as $$\frac{\delta W_{\it eff}}{\delta J}\Bigg|_{J=0} = \langle \phi\rangle.$$ Suppose i have different fields and different sources, labelled for example as $\phi_M,J^M$, $M=1,\cdots,N$. What are the conditions for the functional derivatives to commute, in other words, does the following hold? $$\frac{\delta W_{\it {eff}}}{\delta J^M(x) \delta J^N(y)}=\frac{\delta W_{\it {eff}}}{\delta J^N(x) \delta J^M(y)}?$$ It is not obvious to me whether this is the case, since we are evaluating the functional derivatives, although it makes intuitive sense to have the property.
Somewhat belated, but it is I think safe to consider the functional derivative as an infinite-dimensional generalization of the gradient operator. So they commute whenever gradients commute. Which, on a flat space, is always, I believe.
More explicitly: if you think of every value of $J^M(x), J^N(x)$ as an element of a vector, then a derivative w.r.t. a $J^M(x)$ for particular $x$ amounts to increasing the value of $J^M$ for that $x$ (stepping in a given direction) and seeing how much $W_{\it eff}$ changes. Taking a second derivative w.r.t. $J^N(y)$ amounts to taking a second step in a different direction. Commuting the derivatives reverses the order of operations, but unless your space is curved you'd expect to end up in the same place.
So, unless your domain of integration/definition of your functions (Banach space?) is for some perverse reason itself curved, they should always commute.