We can calculate degeneracy of the ground space of the Toric code in a few different ways.
Using stabilizer formalism
A simple way to find the degeneracy uses the stabilizer formalism. The relevant result says that a stabilizer group generated by $g$ independent generators fixes states in a $2^{n-g}$-dimensional subspace of the $n$-qubit Hilbert space. See Proposition 10.5 on page 458 in <PERSON> & <PERSON>'s "Quantum Computation and Quantum Information".
The stabilizer group of the Toric code on $2k^2$ qubits is generated by $k^2$ site operators $A_s$ and $k^2$ plaquette operators $B_p$ defined as
$$ A_s = \prod_{j \in s} \sigma^x_j \,\,\,\,\,\, B_p = \prod_{j \in p} \sigma^z_j \tag1 $$
(c.f. equation $(1)$ in <PERSON>'s paper). However, since $(\sigma^z_j)^2 = (\sigma^z_j)^2 = I$ we see that any one of the site operators is equal to the product of all the others and similarly for the plaquette operators. Therefore, we have $2k^2-2$ independent generators. By the result cited above, the stabilizer group of the Toric code fixes states in a subspace which has $2^{2k^2 - (2k^2-2)} = 4$ dimensions, i.e. is $4$-fold degenerate.
Using homology
An intuitive way of calculating the dimension of the ground space of the Toric code uses the concept of homology.
The Hamiltonian of the Toric code is defined as
$$ H = -\sum_s A_s - \sum_p B_p $$
where $s$ ranges over sites and $p$ ranges over plaquettes of the lattice (c.f. equation $(4)$ in <PERSON>'s paper). 578
The operators $A_s$ and $B_p$ commute for all $s$ and $p$, so there is a basis that consists of their shared eigenvectors and $H$ is diagonal in this basis. The eigenvalues of $A_s$ and $B_p$ are $+1$ and $-1$, so we see that the lowest eigenvalue of $H$ corresponds to the eigenvectors $|\xi_i\rangle$ of $A_s$ and $B_p$ associated with eigenvalue $+1$. In other words, the ground space of $H$ is
$$ \mathcal{L} = {|\xi_i\rangle\,:\,A_s|\xi_i\rangle = |\xi_i\rangle, B_p|\xi_i\rangle = |\xi_i\rangle\,\text{for all}\,p,s} $$
(c.f. equation $(2)$ in the paper).
Note that this analysis also reveals the gap between closest eigenvalues of $H$. The gap is due to the sign change in the eigenvalue of a site or plaquette operator. However, since sites can only change sign in pairs and plaquettes can only change sign in pairs, we see that the gap is equal to $4$.
Now, for a set $a$ of links of the lattice define
$$ A_a = \prod_{j \in a} \sigma^x_j $$
and similarly for a set $b$ of links of the dual lattice define
$$ B_b = \prod_{j \in b} \sigma^z_j. $$
If $a$ is a closed loop, then $A_a$ commutes with all site and plaquette operators $(1)$. However, only the closed loops which are boundaries of two-dimensional regions can be generated as products of site and plaquette operators. This allows us to define operators which preserve $\mathcal{L}$ while acting non-trivially on it.
The first homology group of the torus has rank $2$ reflecting the fact that there are two types of loops that are not boundaries. Denote the two homologically non-trivial loops of the lattice $\alpha$ (running vertically) and $\beta$ (running horizontally) and similarly for the dual lattice $\alpha'$ (running vertically) and $\beta'$ (running horizontally). 66