Show that $\large X^{q^n}−X \in \mathbb{F}_q[X]$ (with $q = p^k$ for some prime $p \in \mathbb{N}^+$ and $k,n \in \mathbb{N}^+$) is the product of all monic irreducible polynomials $\pi \in \mathbb{F}_q[X]$ with $\deg(\pi) ~|~ n$:
Lemma 1:
$\forall q, n, d \in \mathbb{N}^+: \large q^n \bmod \left(q^d-1\right) = q^{n ~\bmod~ d}$
as $q^d = 1 \bmod \left(q^d-1\right)$
Lemma 2:
$\large \gcd\left(X^{q^n} - X, X^{q^d} - X\right) = X^{q^{\gcd(n, d)}} - X$
(in any polynomial ring over a field, especially in $\mathbb{F}_q[X]$)
For $n = d$ this is obvious, assume w.l.o.g. $n > d$. For all $1 \leq k \in \mathbb{N}$ with $q^n - k(q^d-1) > 0$ (required so all exponents are $\geq 0$):
$$\large X^{q^n} - X = \left(\sum\limits_{i=1}^k X^{q^n -i(q^d-1) - 1}\right)\cdot \left(X^{q^d}-X\right) + \left(X^{q^n -k(q^d-1)} - X\right)$$
As $q^n \bmod \left(q^d-1\right) = q^{n ~\bmod~ d} \neq 0$ ($\exists k: q^n \bmod \left(q^d-1\right) = q^n - k(q^d-1) > 0$):
$$ \large \Rightarrow \left(\large X^{q^n} - X\right) \bmod \left(X^{q^d}-X\right) = \left(X^{q^{n ~\bmod~ d}} - X\right)$$
$$ \large \Rightarrow \gcd\left(X^{q^n} - X, X^{q^d} - X\right) = \gcd\left(X^{q^d} - X, X^{q^{n ~\bmod~ d}} - X\right)$$
I.e. the $\gcd$ modulo reduction is done in the $q$ and $d$ exponents.
Step 1: (Similar to the answer by Thomas Andrews)
Let $\pi$ be a monic irreducible polynomial in $\mathbb{F}_q[X]$ with degree $d = \deg(\pi)$, and $F_{\pi} := \mathbb{F}_q[X]/\left<\pi\right>$ with $\varphi: \mathbb{F}_q[X] \to F_{\pi}$. Show $d ~|~ n \Rightarrow \pi ~|~ \left(X^{q^n}-X\right)$.
As the size of the multiplicative subgroup $\large \left|F_{\pi}^\ast\right| = q^d-1$, it follows that $\large\forall y \in F_{\pi}: y^{q^d-1} = 1, y^{q^d} - y = 0$. $\large y^{q^d} - y = 0$ is also true for $\large y = 0$, i.e. $\large\forall y \in F_{\pi}$.
Therefore $\large \varphi(X^{q^d}-X) = 0 \Rightarrow \exists k \in \mathbb{F}_q[X]: \left(X^{q^d}-X\right) = 0 + k \cdot \pi \Rightarrow \pi ~|~ \left(X^{q^d}-X\right)$
If $d ~|~ n \Rightarrow \large \gcd\left(X^{q^n} - X, X^{q^d} - X\right) = X^{q^d} - X \Rightarrow \pi ~|~ \left(X^{q^n} - X\right)$
Step 2:
$\large f = X^{q^n} - X \in \mathbb{F}_q[X]$ is square-free, as $\large f' = q^n \cdot X^{q^n-1} - 1 = -1$ ($q = 0$ in $\mathbb{F}_q$!), and $\gcd(f, f') = 1$.
If $\exists a, b \in \mathbb{F}_q[X]: f = (a \cdot a) \cdot b $
$\Rightarrow f' = (a\cdot a' + a' \cdot a) \cdot b + (a \cdot a) \cdot b' = a \cdot (a'\cdot b + a'\cdot b + a \cdot b')$
$\Rightarrow a ~|~ \gcd(f, f')$
As $\gcd(f, f') = 1$ there is no $a \in \mathbb{F}_q[X]$ with $\deg(a) \geq 1$ and $a ~|~ \gcd(f, f')$, and $f$ is square-free.
Step 3:
Induction over $n \geq 1$: show that $ \large p_n := X^{q^n} - X \in \mathbb{F}_q[X]$ is the product of all monic irreducible polynomials $\pi \in \mathbb{F}_q[X]$ with $\deg(\pi) ~|~ n$.
We already know that all such $\pi$ are factors of $p_n$ and $p_n$ is square free. Now show that all factors have the required form.
Let $\pi \in \mathbb{F}_q[X]$ be a irreducible polynomial with $d = \deg(\pi)$ and $\pi ~|~ p_n$. Then $\pi ~|~ \gcd(p_n, p_d) = p_{\gcd(n, d)}$.
If $\gcd(n, d) < d$ then (by induction) $\pi \nmid p_{\gcd(n, d)}$ as $d \nmid n$.
Otherwise $\gcd(n, d) = d \Rightarrow d ~|~ n$.