Can the effective vertex for $\gamma\to3\pi$ be derived directly from the anomaly? My question is whether the effective vertex for $\gamma\to3\pi$ can be derived directly from the anomaly (given in the first equation below), in analogy with the $\pi^0\to2\gamma$ vertex? As far as I understand, based on the anomaly one can derive the Wess-Zumino-Witten (WZW) action which contains many more various vertices aside from $\pi^0\to2\gamma$ and $\gamma\to3\pi$. The derivation is however quite tedious and I do not understand it well. On the other hand, the $\pi^0\to2\gamma$ term can be deduced directly from the anomaly, as is shown below. I do not see how to generalize this reasoning to include $\gamma\to3\pi$ vertex but I believe this should be possible. The question is how?
The amplitude for the $\pi^0\to2\gamma$ decay gets contribution from the anomaly in the chiral current
$$\partial_\mu J^{\mu3}_A=-\frac{e^2}{32\pi^2}\tilde{F}^{\mu\nu}F_{\mu\nu}$$
According to my understanding, this can be done as follows. From the effective pion Lagrangian
$$\mathcal{L}=\frac{f_\pi^2}{4}Tr\left(\partial_\mu U \partial^\mu U^{-1}\right)+O(f_\pi^{0})$$
one finds that in terms of the pion fields the chiral current is given as
$$J^{\mu3}_A=f_\pi\partial^\mu\pi^0+O(f_\pi^{-1})$$
Hence, combining this with the anomaly equation one obtains
$$\partial_\mu\partial^\mu\pi^0=\frac{e^2}{32\pi^2f_\pi}\tilde{F}^{\mu\nu}F_{\mu\nu}+O(f_\pi^{-2})$$
This equation of motion would follow from the corresponding term in the effective Lagrangian at order $O(f_{\pi}^{-1})$
$$\Delta \mathcal{L}=\pi^0\frac{e^2}{32\pi^2f_\pi}\tilde{F}^{\mu\nu}F_{\mu\nu}$$
This coupling term indeed produces the correct amplitude for $\pi^0\to2\gamma$.
Alternatively, this term could be derived by expanding the WZW action, which encompasses the effects of anomaly to all orders in $f_{\pi}$. I do not quite understand how the WZW action is derived, and that may be the source of my confusion. Anyway, the WZW action includes much more vertices including for example $\gamma\to 3\pi$
$$\Delta\mathcal{L}\propto \epsilon^{\mu\alpha\beta\gamma}\epsilon^{abc}A_\mu \partial_\alpha\pi^a\partial_\beta\pi^b\partial_\gamma\pi^c$$
Contrary to the vertex for $\pi\to2\gamma$ this one contains a single photon. It is not clear to me how this term can be derived from the anomaly. However, I believe that this should be possible. Is it?

A side question: $\pi^0\to2\gamma$ piece is clearly not isospin-invariant, since it contains only the $\pi^0$  field. Can one think of the WZW action as the one restoring the symmetry? Can one reconstruct the full WZW action from this single piece?
Edit: I've run into confusion here. As the axial anomaly shows, the isospin symmetry is broken by the electromagnetic interactions, so my suggestion above makes no sense.
 A: Motivation for the Wess-Zumino term
There are many similar terms with mesons which capture anomalous chiral structure of underlying QCD. How to capture all of them? The answer is given by the theorem, which states following: in general, non-abelian anomaly in 4 dimensions is related Chern-Simons character in 5 dimensions, which is called the Wess-Zumino term. Such theorem also is an elegant way to derive the Wess-Zumino term.
So instead of looking for all terms like $\pi^{0}F_{EM}\tilde{F}_{EM}$ by direct application of anomaly equations it's preferable to gauge the Wess-Zumino term, since it is straightforward. Many anomalous terms arise not from triangle diagram, but, for example, from anomalous pentagonal diagram, which leads to great complication of calculations based on modification of naive Ward identities for given class of diagrams.
How to gauge the Wess-Zumino term? There is an old simple way which is called trial and errors method. Namely, let's start from the Wess-Zumino free term:
$$
N_{c}\Gamma_{WZ} = i\frac{N_{c}}{240 \pi^{2}}\int d^{5}x\epsilon^{ijklm}\text{Tr}\left[L_{i}L_{j}L_{k}L_{l}L_{m}\right],
$$
where $L_{i} \equiv U\partial_{i}U^{-1}$ and $U$ is the goldstone bosons matrix.
Since $U$ here acts in 5-dimensional manifold, there is no direct way to gauge the Wess-Zumino term by simple elongation of derivative. However, we know the explicit form of $U_{EM}(1)$ gauge variation of $U$ field, namely
$$
U(x) \to e^{iQ\epsilon(x)}Ue^{-iQ\epsilon (x)} \Rightarrow \delta U = i\epsilon (x)[Q, U]
$$
Then
$$
\delta L_{i} = i\epsilon (x)[Q, L_{i}] + i(\partial_{i}\epsilon (x))UQU^{-1} - i\partial_{i}\epsilon (x)Q,
$$
and
$$
N_{c}\Delta \Gamma_{WZ} = -\frac{iN_{c}}{48 \pi^{2}}\int d^{5}x\partial_{i}\left[ \partial_{i}\epsilon(x)\epsilon^{ijklm}\text{Tr}\left[ Q(T_{k}T_{l}T_{m} + L_{k}L_{l}L_{m})\right]\right] =
$$
$$
=-\frac{i}{48 \pi^2}\int d^{4}x\epsilon^{\mu \nu \alpha \beta}\partial_{\mu}\epsilon (x) \text{Tr}\left[ QT_{\nu}T_{\alpha}T_{\beta} + L_{\nu}L_{\alpha}L_{\beta}\right] \equiv -\int d^{4}x\partial_{\mu}\epsilon(x)J^{\mu},
$$
where $T_{i} \equiv (\partial_{i}U^{-1})U$ and 
$$
J^{\mu} \equiv \frac{iN_{c}}{48 \pi^{2}}\epsilon^{\mu \nu \alpha \beta}\text{Tr}\left[ Q(T_{\nu}T_{\alpha}T_{\beta} + L_{\nu}L_{\alpha}L_{\beta})\right]
$$
For first step, we may modify the Wess-Zumino term by adding
$$
-\int d^{4}xA_{\mu}J^{\mu},
$$
where $A_{\mu}$ is EM 4-potential.
Since $\delta J_{\mu}$ isn't zero, we need to add also the part whose variation coincides term with $\delta J_{\mu}$. Thus we obtain that the Gauged Wess-Zumino term is
$$
\tag 1 \tilde{\Gamma}_{WZ} = N_{c}\Gamma_{WZ} - \int d^{4}xA_{\mu}J^{\mu} + 
$$
$$
+\frac{iN_{c}}{24 \pi^{2}}\int d^{4}x\epsilon^{\mu \nu \alpha \beta}(\partial_{\mu}A_{\nu})A_{\alpha}\text{Tr}\left[ Q^{2}(\partial_{\beta}U)U^{-1} + Q^2Q^{-1}\partial_{\beta}U + QUQU^{-1}(\partial_{\beta}U)U^{-1}\right]
$$
For QCD with $u, d$ quarks, $Q = \text{diag}\left(\frac{2}{3}, -\frac{1}{3}, \right)$. By using 
$$
U = e^{i\frac{\sqrt{2}\pi_{a}t_{a}}{f_{\pi}}}, \quad \pi_{a}t_{a} = \begin{pmatrix} \frac{\pi^{0}}{\sqrt{2}} & \pi^{-} \\ \pi^{+} & \frac{\pi^{0}}{\sqrt{2}}\end{pmatrix},
$$
we may get $\pi^{0} \to 2\gamma$ vertex from the last summand of $(1)$ and $\gamma \to 3 \pi$ vertex from the second one. Of course, by enlarging the number of quarks to 3, we get more anomalous terms.
