There is a subtle but heavy abuse of the notation that renders many of the steps confusing. Let's address this issue by going back to the definitions of matrix multiplication, transposition, traces, and derivatives. For those wishing to omit the explanations, just jump to the last section "Putting It All Together" to see how short and simple a rigorous demonstration can be.
Notation and Concepts
Dimensions
For the expression $ABA^\prime C$ to make sense when $A$ is an $m\times n$ matrix, $B$ must be a (square) $n\times n$ matrix and $C$ must be an $m\times p$ matrix, whence the product is an $m\times p$ matrix. In order to take the trace (which is the sum of diagonal elements, $\operatorname{Tr}(X)=\sum_i X_{ii}$), then $p=m$, making $C$ a square matrix.
Derivatives
The notation "$\nabla_A$" appears to refer to the derivative of an expression with respect to $A$. Ordinarily, differentiation is an operation performed on functions $f:\mathbb{R}^N\to\mathbb{R}^M$. The derivative at a point $x\in \mathbb{R}^N$ is a linear transformation $Df(x):\mathbb{R}^N\to\mathbb{R}^M$. Upon choosing bases for these vector spaces, such a transformation can be represented as an $M\times N$ matrix. That is not the case here!
Matrices as vectors
Instead, $A$ is being considered as an element of $\mathbb{R}^{mn}$: its coefficients are being unrolled (usually either row by row or column by column) into a vector of length $N=mn$. The function $f(A)=\operatorname{Tr}(ABA^\prime C)$ has real values, whence $M=1$. Consequently, $Df(x)$ must be a $1\times mn$ matrix: it's a row vector representing a linear form on $\mathbb{R}^{mn}$. Howver, the calculations in the question use a different way of representing linear forms: their coefficients are rolled back up into $m\times n$ matrices.
The trace as a linear form
Let $\omega$ be a constant $m\times n$ matrix. Then, by definition of the trace and of matrix multiplication,
$$\eqalign{
\operatorname{Tr}(A\omega^\prime) &= \sum_{i=1}^m(A\omega^\prime)_{ii} = \sum_{i=1}^m\left(\sum_{j=1}^n A_{ij}(\omega^\prime)_{ji}\right) = \sum_{i,j} \omega_{ij}A_{ij}
}$$
This expresses the most general possible linear combination of the coefficients of $A$: $\omega$ is a matrix of the same shape as $A$ and its coefficient in row $i$ and column $j$ is the coefficient of $A_{ij}$ in the linear combination. Because $\omega_{ij}A_{ij}=A_{ij}\omega_{ij}$, the roles of $\omega$ and $A$ may switched, giving the equivalent expression
$$\sum_{i,j} \omega_{ij}A_{ij} = \operatorname{Tr}(A\omega^\prime) = \operatorname{Tr}(\omega A^\prime).\tag{1}$$
By identifying a constant matrix $\omega$ with either of the functions $A\to \operatorname{Tr}(A \omega^\prime)$ or $A\to \operatorname{Tr}(\omega A^\prime)$, we may represent linear forms on the space of $m\times n$ matrices as $m\times n$ matrices. (Do not confuse these with derivatives of functions from $\mathbb{R}^n$ to $\mathbb{R}^m$!)
Computing a Derivative
The definition
Derivatives of many of the matrix functions encountered in statistics are most easily and reliably computed from the definition: you don't really need to resort to complicated rules of matrix differentiation. This definition says that $f$ is differentiable at $x$ if and only if there is a linear transformation $L$ such that
$$f(x+h) - f(x) = Lh + o(|h|)$$
for arbitrarily small displacements $h\in \mathbb{R}^N$. The little-oh notation means that the error made in approximating the difference $f(x+h)-f(x)$ by $Lh$ is arbitrarily smaller than the size of $h$ for sufficiently small $h$. In particular, we may always ignore errors that are proportional to $|h|^2$.
The calculation
Let's apply the definition to the function in question. Multiplying, expanding, and ignoring the term with a product of two $h$'s in it,
$$\eqalign{
f(A+h)-f(A) &= \operatorname{Tr}((A+h)B(A+h)^\prime C) - \operatorname{Tr}(ABA^\prime C) \\
&= \operatorname{Tr}(hBA^\prime C) +\operatorname{Tr}(ABh^\prime C) + o(|h|).\tag{2}
}$$
To identify the derivative $L=Df(A)$, we must get this into the form $(1)$. The first term on the right is already in this form, with $\omega = BA^\prime C$. The other term on the right has the form $\operatorname{Tr}(Xh^\prime C)$ for $X=AB$. Let's write this out:
$$\operatorname{Tr}(Xh^\prime C) = \sum_{i=1}^m\sum_{j=1}^n\sum_{k=1}^m X_{ij} h_{kj} C_{ki} = \sum_{i,j,k}h_{kj} \left(C_{ki}X_{ij}\right) =\operatorname{Tr}((CX)h^\prime).\tag{3}$$
Recalling $X=AB$, $(2)$ can be rewritten
$$f(A+h) - f(A) = \operatorname{Tr}(h\, BA^\prime C\,) + \operatorname{Tr}(CAB\, h^\prime\,)+o(|h|).$$
It is in this sense that we may consider the derivative of $f$ at $A$ to be $$Df(A) = (BA^\prime C)^\prime + CAB = C^\prime A B^\prime + CAB,$$ because these matrices play the roles of $\omega$ in the trace formulas $(1)$.
Putting It All Together
Here, then, is a complete solution.
Let $A$ be an $m\times n$ matrix, $B$ an $n\times n$ matrix, and $C$ an $m\times m$ matrix. Let $f(A) = \operatorname{Tr}(ABA^\prime C)$. Let $h$ be an $m\times n$ matrix with arbitrarily small coefficients. Because (by identity $(3)$) $$\eqalign{f(A+h) - f(A) &= \operatorname{Tr}(hBA^\prime C) +\operatorname{Tr}(ABh^\prime C) + o(|h|) \\
&=\operatorname{Tr}(h(C^\prime A B^\prime)^\prime + (CAB)h^\prime) + o(|h|),}$$ $f$ is differentiable and its derivative is the linear form determined by the matrix $$C^\prime A B^\prime + CAB.$$
Because this takes only about half the work and involves only the most basic manipulations of matrices and traces (multiplication and transposition), it has to be considered a simpler--and arguably more perspicuous--demonstration of the result. If you really want to understand the individual steps in the original demonstration, you might find it fruitful to compare them to the calculations shown here.
