In this answer we shall use the Schwinger action principle to prove the CCRs, as OP requested. Let us for simplicity consider the Hamiltonian formulation of bosonic point mechanics. (It is in principle possible to generalize to the Lagrangian formulation, field theory, & fermionic variables.)
We start with the Hamiltonian action
$$\begin{align} S_H(t_i,t_f)~=~&\int_{t_i}^{t_f}\!dt~ L_H, \cr L_H~=~&\sum_{k=1}^np_k \dot{q}^k - H . \end{align} \tag{1}$$
The Hamiltonian phase space path integral reads$^1$
$$\begin{align} \langle q_f,t_f | q_i, t_i \rangle
~=~&\langle q_f,t |\exp\left\{-\frac{i}{\hbar} \hat{H} \Delta t \right\}| q_i, t \rangle \cr
~=~&\int_{q(t_i)=q_i}^{q(t_f)=q_f} \!{\cal D}q~{\cal D}p~\exp\left\{\frac{i}{\hbar} S_H(t_i,t_f)\right\}, \end{align} \tag{2} $$
$$ \langle q_f,t_f | \hat{F}|q_i, t_i \rangle
~=~ \int_{q(t_i)=q_i}^{q(t_f)=q_f} \!{\cal D}q~{\cal D}p~F~\exp\left\{\frac{i}{\hbar} S_H(t_i,t_f)\right\} . \tag{3} $$
If we change infinitesimally the action (1), we derive the Schwinger action principle:
$$\begin{align} \frac{\hbar}{i}\delta\langle q_f,t_f &| q_i, t_i \rangle \cr
~\stackrel{(2)}{=}~& \int_{q(t_i)=q_i}^{q(t_f)=q_f} \!{\cal D}q~{\cal D}p~ \delta_0 S_H(t_i,t_f)~\exp\left\{\frac{i}{\hbar} S_H(t_i,t_f)\right\}\cr
~\stackrel{(3)}{=}~&\langle q_f,t_f | \widehat{\delta_0 S}_H(t_i,t_f)|q_i, t_i \rangle . \end{align} \tag{4} $$
Similarly, we are interested in calculating the change to
$$\begin{align} \langle q^{\prime}_i,t_i |& \hat{F}(t_f)|q^{\prime\prime}_i, t_i \rangle \cr
~=~& \int \! dq^{\prime}_f~dq^{\prime\prime}_f~
\langle q^{\prime}_i,t_i |q^{\prime}_f,t_f \rangle ~
\langle q^{\prime}_f,t_f | \hat{F}(t_f)|q^{\prime\prime}_f, t_f \rangle~\langle q^{\prime\prime}_f,t_f |q^{\prime\prime}_i,t_i \rangle ,\end{align} \tag{5} $$
where the initial states $|q, t_i \rangle$ and the middle bracket on the rhs. of eq. (5) are unaffected by the (later) change of action in the (later) time interval $[t_i,t_f]$.
We calculate the change
$$\begin{align} \langle q^{\prime}_i,t_i | &\delta\hat{F}(t_f)|q^{\prime\prime}_i,t_i\rangle\cr
~=~& \delta\langle q^{\prime}_i,t_i | \hat{F}(t_f)|q^{\prime\prime}_i,t_i\rangle\cr
~\stackrel{(5)}{=}~&
\int \! dq^{\prime}_f~dq^{\prime\prime}_f~
\delta\langle q^{\prime}_i,t_i |q^{\prime}_f,t_f \rangle ~
\langle q^{\prime}_f,t_f | \hat{F}(t_f)|q^{\prime\prime}_f, t_f \rangle~\langle q^{\prime\prime}_f,t_f |q^{\prime\prime}_i,t_i \rangle\cr
~+~&\int \! dq^{\prime}_f~dq^{\prime\prime}_f~
\langle q^{\prime}_i,t_i |q^{\prime}_f,t_f \rangle ~
\langle q^{\prime}_f,t_f | \hat{F}(t_f)|q^{\prime\prime}_f, t_f \rangle~\delta\langle q^{\prime\prime}_f,t_f |q^{\prime\prime}_i,t_i \rangle\cr
~\stackrel{(4)}{=}~&\frac{i}{\hbar}\langle q^{\prime}_i,t_i | [\hat{F}(t_f),\widehat{\delta_0 S}_H(t_i,t_f)]|q^{\prime\prime}_i,t_i\rangle ,\end{align} \tag{6}$$
or equivalently, as an operator identity
$$ \delta\hat{F}(t_f)~\stackrel{(6)}{=}~\frac{i}{\hbar} [\hat{F}(t_f),\widehat{\delta_0 S}_H(t_i,t_f)]. \tag{7}$$
In other words, the cause $\delta_0$ leads to the effect $\delta$.
We next go to the adiabatic limit $\Delta t:=t_f-t_i\to 0$ where the symplectic term of the action (1) dominates over the Hamiltonian term, i.e. we can effectively remove the Hamiltonian $H$ from the calculation, i.e. operators in the Heisenberg picture can be treated as time-independent in this limit.$^2$
$$\begin{align} \frac{\hbar}{i} \delta\hat{q}^k(t_f)
~\stackrel{(7)}{\simeq}~&
[\hat{q}^k(t_f),\widehat{\delta_0 S}_H(t_i,t_f)]\cr
~\stackrel{(1)}{\simeq}~&
\left[\hat{q}^k(t_f), \sum_{\ell=1}^n \hat{p}_{\ell}(t_f)~
\delta_0\{ \hat{q}^{\ell}(t_f) - \hat{q}^{\ell}(t_i)\} \right]\cr
~=~&\sum_{\ell=1}^n [\hat{q}^k(t_f), \hat{p}_{\ell}(t_f)]~\delta_0\hat{q}^{\ell}(t_f) . \end{align} \tag{8}$$
Since there is effectively no Hamiltonian, the cause & effect must cancel:
$$\delta_0\hat{q}^{\ell}(t_f) + \delta\hat{q}^{\ell}(t_f)~=~0. \tag{9}$$
Eq. (8) is only possible if we have the equal-time CCR
$$\begin{align} [\hat{q}^k(t), \hat{p}_{\ell}(t)]~=~&i\hbar \delta_{\ell}^k~{\bf 1} , \cr k,\ell~\in~&\{1,\ldots, n\},\end{align} \tag{10}$$
as OP wanted to show. $\Box$
--
$^1$ Here we will work in the Heisenberg picture with time-independent ket and bra states, and time-dependent operators
$$\begin{align} \hat{F}(t_f)~=~&\exp\left\{\frac{i}{\hbar} \hat{H} \Delta t \right\}
\hat{F}(t_f)\exp\left\{-\frac{i}{\hbar} \hat{H} \Delta t \right\}, \cr \Delta t~:=~&t_f-t_i. \end{align} \tag{11}$$
Moreover, $| q, t \rangle $ are instantaneous position eigenstates in the Heisenberg picture,
$$\begin{align} \hat{q}^k(t)| q, t \rangle~=~&q^k| q, t \rangle, \cr k~\in~&\{1,\ldots, n\}.\end{align} \tag{12}$$
see e.g. J.J. Sakurai, Modern Quantum Mechanics, Section 2.5. Such states (12) are only well-defined for commuting observables
$$ [\hat{q}^k(t), \hat{q}^{\ell}(t)]~=~0 , \qquad k,\ell~\in~\{1,\ldots, n\}, \tag{13}$$
so we are not going to derive the CCR (13). Rather (13) is an assumption with this proof.
Finally for completeness, let us mention that instead of instantaneous position eigenstates $| q, t \rangle $, we could use instantaneous momentum eigenstates $| p, t \rangle $, but again we would have to assume the corresponding CCR
$$ [\hat{p}_k(t), \hat{p}_{\ell}(t)]~=~0 , \qquad k,\ell~\in~\{1,\ldots, n\}. \tag{14}$$
Assumptions (13) & (14) can be avoided by using other methods. E.g. the CCRs (10), (13) & (14) also follows from the Peierls bracket.$^3$
$^2$ Notation:
The $\approx$ symbol means equality modulo eqs. of motion.
The $\sim$ symbol means equality modulo boundary terms.
The $\simeq$ symbol means equality in the adiabatic limit, where the Hamiltonian $H$ can be ignored.
$^3$ The correspondence principle between quantum mechanics and classical mechanics states that the equal-time CCR
$$\begin{align} [\hat{z}^I(t), \hat{z}^K(t)]
~\stackrel{(16)}{=}~&i\hbar\omega^{IK}~{\bf 1} , \cr I,K~\in~&\{1,\ldots, 2n\}, \end{align} \tag{15}$$
is $i\hbar$ times the equal-time Peierls bracket
$$\begin{align} \{z^I(t), z^K(t^{\prime})\}
~\stackrel{(17)}{\simeq}~&\omega^{IK} , \cr I,K~\in~&\{1,\ldots, 2n\}. \end{align} \tag{16}$$
The Peierls bracket is defined as
$$\begin{align} \{ F,G \}~:=~&\iint_{[t_i,t_f]^2}\!dt~dt^{\prime}~\sum_{I,K=1}^{2n} \frac{\delta F }{\delta z^I(t)}~G^{IK}_{\rm ret}(t,t^{\prime})~\frac{\delta G }{\delta z^K(t^{\prime})} \cr
&- (F\leftrightarrow G)\cr
~\stackrel{(22)}{\simeq}~&\iint_{[t_i,t_f]^2}\!dt~dt^{\prime}~\sum_{I,K=1}^{2n} \frac{\delta F }{\delta z^I(t)}~\omega^{IK}~\frac{\delta G }{\delta z^K(t^{\prime})}. \end{align} \tag{17}$$
Classically, if the Hamiltonian action (1) is recast in symplectic notation
$$\begin{align} S_H(t_i,t_f)~\sim~&\int_{[t_i,t_f]}\!dt\left(\frac{1}{2}\sum_{I,K=1}^{2n} z^I ~\omega_{IK}~ \dot{z}^K-H\right)\cr
~\sim~&\frac{1}{2} \iint_{[t_i,t_f]^2}\!dt~dt^{\prime}\sum_{I,K=1}^{2n} z^I(t)~\omega_{IK}~ \delta^{\prime}(t\!-\!t^{\prime}) ~z^K(t^{\prime})\cr &-\int_{[t_i,t_f]}\!dt~H, \end{align} \tag{18}$$
and changed infinitesimally, then the classical solution $z^I(t)$ is also changed infinitesimally $\delta z^I(t)$, to ensure that the deformed EL eqs.
$$0~\approx~ \frac{\delta}{\delta z^I(t)} (S_H+\delta_0 S_H)[z+\delta z] \tag{19}$$
are satisfied. (We are ignoring boundary terms everywhere in this calculation.) In other words
$$ \frac{\delta ~\delta_0 S_H(t_i,t_f)}{\delta z^I(t)}~\approx~-\int_{[t_i,t_f]}\! dt^{\prime} \sum_{K=1}^{2n}H_{IK}(t,t^{\prime})~\delta z^K(t^{\prime}), \tag{20}$$
where the Hessian is
$$\begin{align} H_{IK}(t,t^{\prime})
~:=~& \frac{\delta^2 S_H(t_i,t_f)}{\delta z^I(t) \delta z^K(t^{\prime})}\cr
~=~&\omega_{IK}~\delta^{\prime}(t\!-\!t^{\prime}) -\partial_I\partial_K H ~\delta(t\!-\!t^{\prime})\cr
~\simeq~&\omega_{IK}~\delta^{\prime}(t\!-\!t^{\prime}). \end{align} \tag{21}$$
Then the retarded Green's function simplifies to
$$G^{IK}_{\rm ret}(t,t^{\prime})
~\stackrel{(21)+(23)}{\simeq}~\omega^{IK}~\theta(t\!-\!t^{\prime}),\tag{22}$$
$$
\int_{[t_i,t_f]}\! dt^{\prime} \sum_{J=1}^{2n}H_{IJ}(t,t^{\prime})~G^{JK}_{\rm ret}(t^{\prime},t^{\prime\prime})
~=~\delta_I^K~\delta^{\prime}(t\!-\!t^{\prime\prime}). \tag{23}$$
Therefore we derive the classical analogue of the operator identity (7)
$$ \delta F(t_f)~\stackrel{(25)}{=}~\{\delta_0 S_H(t_i,t_f), F(t_f)\} \tag{24} $$
from
$$\begin{align} \delta z^I(t_f)
~\stackrel{(20)+(23)}{=}&-\int_{[t_i,t_f]}\! dt^{\prime} \sum_{K=1}^{2n} G^{IK}_{\rm ret}(t_f,t^{\prime})~\frac{\delta ~\delta_0 S_H(t_i,t_f)}{\delta z^K(t^{\prime})}\cr
~\stackrel{(17)}{=}~~&\{\delta_0 S_H(t_i,t_f), z^I(t_f)\} \cr
~\simeq~~~&\left\{ \sum_{J=1}^{2n} z^J(t_f)~\omega_{JK} \delta_0z^K(t_f) , z^I(t_f)\right\}\cr
~=~~~&-\sum_{J=1}^{2n} \delta_0 z^J(t_f)~\omega_{JK} \{z^K(t_f) , z^I(t_f)\}\cr
~\stackrel{(16)}{\simeq}~~&-\delta_0z^I(t_f),\end{align} \tag{25}$$
which in turn is consistent with the fact that cause & effect must cancel:
$$\delta_0z^I(t_f) + \delta z^I(t_f)~=~0, \tag{26}$$
since there is effectively no Hamiltonian.