There are two ways to quantize a gauge theory.
1) First quantize all degrees of freedom and then reduce the gauge freedom by imposing conditions on the quantum Hilbert space.
2) First reduce the gauge redundancy then quantize the reduced phase space.
For realistic gauge theories in $4D$, there are problems in implementing either of these methods. The first method would require gauge fixing and the Faddeev-Popov procedure, but actually there is no good gauge fixing due to the Gribov ambiguity. The second method would result a very complicated non-flat infinite dimensional phase space with singularities.
In the case of the pure $3D$ Chern-Simons however, the second method results a finite dimensional quantizable phase space. In the Abelian (and also Non-Abelian) Chern-Simons theory on the torus, this construction can be performed explicitly.
Recall that the phase space of a quantum theory can be identified with the space of its classical solutions. In the case of the Chern Simons theory the classical solution satisfy:
$$\mathbf{F}=0$$
($\mathbf{F}$ is the gauge field strength). In flat space, this condition would mean that the solution is a pure gauge and the removal of the gauge redundancy would mean that the phase space is one point. However on the torus there is another solution which is not a pure gauge: The constant solution, because on one hand it satisfies the classical equation of motion and on the other hand it is normalizable. Thus a classical solution on the torus would have the general form:
$$\mathbf{A(x)}=\frac{\mathbf{\alpha}}{2\pi}+\mathbf{\nabla}{\phi(\mathbf{x})}$$
Where $\mathbf{\alpha}$ is a constant.
Remembering that the Poisson brackets obtained from the classical Abelian Chern Simons Lagrangian read:
$$\{A_i, A_j \} = \frac{2 \pi i \epsilon_{ij}}{k} \delta^2(\mathbf{x}-\mathbf{y})$$
And observing that
$$\mathbf{\alpha} = \oint\mathbf{A}$$
(The integration is over a single turn). Then by integrating the Poisson bracket equation over its two circle generators, we obtain the following Poisson brackets for $\alpha$
$$\{\alpha_i, \alpha_j \} = \frac{2 \pi i \epsilon_{ij}}{k} $$
$\alpha$ defines a two dimensional phase space. Quantizing this space would mean to find functions on the two dimensional phase space satisfying these Poisson brackets. In this case this can be accomplished by inspection:
$$ \alpha_i = p_i - \frac{2 \pi i \epsilon_{ij}}{k} q_i$$
This is just the momentum operator in the presence of a constant magnetic field $$\mathbf{B} = \frac{2 \pi }{k} \mathbf{\hat{z}}$$ in the $z$-direction. ("Magnetic translation operator"). Thus the quantization of the Wilson loop along an integral number of loops:
$$ W(\mathbf{l}) = e^{i \oint_{\mathbf{l}}\mathbf{A}} = e^{i \mathbf{\alpha}. \mathbf{l}}$$
(i.e., this time we do not restrict the integration to a single turn:) $\mathbf{l} = 2 \pi (m, n)$ where $m$ and $n$ are integers.
It is not hard using the canonical commutation relations (Weyl algebra) to obtain:
$$ W(\mathbf{l}_1) W(\mathbf{l}_2) = e^{i \mathbf{B} .( \mathbf{l}_1 \times \mathbf{l}_2)} W(\mathbf{l}_2) W(\mathbf{l}_1) $$
In particular:
$$ W(\mathbf{\hat{x}}) W(\mathbf{\hat{y}}) = e^{i \frac{2 \pi}{k}} W(\mathbf{\hat{y}}) W(\mathbf{\hat{x}}) $$
It is worthwhile to mention that the operators $\alpha$ cannot be chosen as observables because they are not gauge invariant. Thus we seek a Hilbert space representation of $W(\mathbf{\hat{x}})$ and $W(\mathbf{\hat{y}})$. One way to find such a representation is to postulate the existence of a state (vacuum) $|0\rangle$ upon which we can generate the whole Hilbert space by the action of the operators $W(\mathbf{\hat{x}})$ and $W(\mathbf{\hat{y}})$.
First, we observe that, without loss of generality, we can choose $W(\mathbf{\hat{x}})$ as diagonal. Next we observe that
$$ W(\mathbf{\hat{x}}) W(\mathbf{\hat{y}})^k = W(\mathbf{\hat{y}})^k W(\mathbf{\hat{x}}) $$
Thus $ W(\mathbf{\hat{y}})^k $ commutes with the diagonal $W(\mathbf{\hat{x}})$, thus they must be diagonal, since it commutes also with $W(\mathbf{\hat{y}})$, therefore with any sequence of products of the two operators. Thus this operator must be represented by the unit operator:
$$ W(\mathbf{\hat{y}})^k = \mathbf{I}$$
Since $ W(\mathbf{\hat{x}})$ is diagonal, it does not alter the vectors it acts on, thus, a basis of the Hilbert space must be given by:
$$\mathrm{Span} \{ |0\rangle$, $ W(\mathbf{\hat{y}}) |0\rangle, W(\mathbf{\hat{y}})^2 |0\rangle,..., W(\mathbf{\hat{y}})^{k-1} |0\rangle \}$$
This implies that the Hilbert space is finite dimensional, and if we denote:
$$ |n\rangle = W(\mathbf{\hat{y}})^n |0\rangle$$
We obtain the torus algebra action on the Hilbert space:
$$ W(\mathbf{\hat{x}}) |n\rangle = e^{i \frac{2 \pi n}{k}} |n\rangle $$
$$ W(\mathbf{\hat{y}}) |n\rangle = |n+1\rangle $$
This is not the only way to obtain this result. This Hilbert space can be obtained more naturally as the space of Jacobi theta functions in the coordinate representation. This would require to use the techniques of Berezin quantization or Coherent state quantization, please see, for example, section 6 of the following article by: Spradlin, and Volovich.