I don't have a reference, but here's the next best thing: a proof. 

First, let's fix some notation. Let $P$ be a p-group, $G$ it's group of automorphisms, and $\Phi(P) = P^p[P,P]$ it's Frattini subgroup. Define inductively $\Phi^k(P)$ as $\Phi(\Phi^{k-1}(P)).$ 

The subgroups $\Phi^k(P)$ form a decreasing chain of characteristic subgroups which exhaust $P.$

Let $$G_k := \ker(G \rightarrow Aut(P/\Phi^k(P)))$$
and
$$G_k' := \ker(G \rightarrow Aut(\Phi^k(P)/\Phi^{k+1}(P))).$$


Then the subgroups $G_k$ form an increasing chain of normal subgroups of $G$ which exhaust $G.$ 

Let $d_k = \dim_{\mathbb{F}_p}(\Phi^k(P)/\Phi^{k+1}(P)).$ The group $P$ can be generated by $d_1$ elements. Choose a generating set $g_1 ... g_{d_1}$ and consider the map from $G_k \cap G_k'/G_{k+1}$ to $(\Phi^k(P)/\Phi^{k+1}(P))^{d_1}$

given by 

$$\sigma \mapsto (\sigma(g_i)g_i^{-1})_{i=1}^{d_1}.$$

This map is an injective group homomorphism.

On the other hand $G_k/(G_{k} \cap G_{k}')$ injects into  $Aut(\Phi^k(P)/\Phi^{k+1}(P)) \cong GL_{d_i}( \mathbb{F}_p).$

Note that if $p^n = |P|$ and $r = max\{d_i : d_i \neq 0\},$ then $\displaystyle\sum_{i=1}^r d_i = n.$ It follows that the order $$|G| = \displaystyle\prod_{k=0}^{r}|G_k/G_{k+1}| = \displaystyle\prod_{k=0}^{r} |(G_k \cap G_k')/G_{k+1}||(G_k \cap G_k')/G_{k+1}|$$ divides 


$$\displaystyle\prod_{s=1}^{d_1-1} (p^{d_1} - p^s)\displaystyle\prod_{k=1}^{r} p^{d_kd_1}\displaystyle\prod_{s=1}^{d_k-1} (p^{d_k} - p^s)$$ 

which divides 

$$\displaystyle\prod_{k=0}^{n-1} (p^n - p^k).$$