With some technical assumption, in fact, CCR arises from requiring compatibility of Lagrangian/Hamiltonian formalism and standard formulation of Quantum Theory. I will not be completely rigorous below and I only indicate the way to follow to achieve the wanted result.
Consider the Hamiltonian of the theory, it reads
$$H(t_0) = \frac{1}{2} \int_{\mathbb{R}^3, t=t_0} \Pi_\mu \Pi^\mu + \vec{\nabla} A_\mu \cdot \vec{\nabla} A_\mu d^3x\tag{1}$$
Form Heisenberg evolution of operators we are committed to assume that
$$\partial_t A_\mu(t,x) = i[H(t), A_\mu(t,x)]\tag{2}\:.$$
On the other hand, from Hamilton equations, or directly form the definition of conjugated momentum in Lagrangian formulation, we have
$$\Pi_\mu(t,x) = \partial_t A_\mu(t,x)\tag{3}\:.$$
So, putting (1) and (2) together, it must be
$$\int [\Pi_\mu(t,x) \Pi^\mu(t,x), A^\nu(t,y) ] d^3x+ \int [\vec{\nabla} A_\mu(t,x) \cdot \vec{\nabla} A_\mu(t,x), A^\nu(t,y)] d^3x = -2i\partial_t A^\nu(t,y) $$
If we now assume that
H1. measurements at different positions at the same fixed time of generally different components of $A_\mu$ are compatible in quantum sense,
i.e.,
$[A_\mu(t,x), A^\nu(t,y)]=0$,
it remains
$$\int [\Pi_\mu(t,x) \Pi^\mu(t,x), A^\nu(t,y) ] d^3x= -2i\partial_t A^\nu(t,y) $$
that is
$$\int \Pi_\mu(t,x)\:[\Pi^\mu(t,x), A^\nu(t,y) ] d^3x
+ \int [\Pi_\mu(t,x), A^\nu(t,y) ]\: \Pi^\mu(t,x) d^3x= -2i\partial_t A^\nu(t,y) $$
If now we further assume that
H2. $[\Pi_\mu(t,x), A^\nu(t,y) ]$ is a number
$[\Pi_\mu(t,x), A^\nu(t,y) ]= c(t,x,y)I$ so that it commutes with operators,
we have
$$\int \Pi_\mu(t,x)c(t,x,y) d^3x
= -i\partial_t A^\nu(t,y)\:. $$ From (3)
$$\int \Pi_\mu(t,x)c(t,x,y) d^3x
= -i\Pi_\nu(t,y)\:. $$
This means
$$\int \Pi_\mu(t,x)\left(c(t,x,y) + i\delta(x,y) \delta^\mu_\nu \right)d^3x =0\:.\tag{5}$$
This identity must be interpreted within the fixed-time smearing procedure (also the previous lines should be interpreted this way, but here I make explicit the formalism since a crucial further hypothesis needs which is stated with this formalism): The operators $A(t,x)$ and $\Pi(t,x)$ have to be smeared with smooth compactly supported functions $f : \mathbb{R}^3 \to \mathbb{R}$ giving rise to the smeared field operators, the ones with mathematical sense.
$$A(t,f) := \int A(t,x) f(x) d^3x\:,\quad \Pi(t,f) := \int \Pi(t,x) f(x) d^3x$$
For instance
$$[A(t,x), \Pi(t,y)] = i \delta(x-y)I$$
has to be interpreted as a short way to write
$$[A(t,f), \Pi(t,g)] = i \int f(x)g(x) d^3x \:I$$
This way (5) actually means that, for every smooth compactly supported function $f : \mathbb{R}^3 \to \mathbb{R}$,
$$\int d^3y \int \Pi_\mu(t,x)\left(c(t,x,y)f(y) + i\delta(x,y) f(y)\delta^\mu_\nu \right)d^3x =0 $$
In other words
$$\Pi_\mu\left(t, \int c(t,\cdot,y)f(y) d^3y +i \delta^\mu_\nu f\right) =0\tag{6}$$
The last hypothesis I demand is that
H3. the operator valued distribution $f \mapsto \Pi_\mu(t,f)$ vanishes if and only if $f=0$.
Assuming this, (6) implies
$$ \int c(t,x,y)f(y) d^3y +i \delta^\mu_\nu f(x) =0$$
for every said function $f$, which means
$$c(t,x,y) = -i\delta^\mu_\nu \delta(x-y)\:,$$
as wanted.
