Starting with some background information from Wikipedia, we have that under time reversal the position is unchanged while the momentum changes sign.
In quantum mechanics we can express the action of time reversal on these operators as $\Theta\,\mathbf{x}\,\Theta^\dagger = \mathbf{x}$ and $\Theta\,\mathbf{p}\,\Theta^\dagger = -\mathbf{p}$. It is worth mentioning here that the time reversal operator, $\Theta$, is anti-unitary, which allows it to be expressed as $\Theta = UK$ where $U$ is unitary and $K$ is the complex conjugation operator.
As for the creation/annihilation operators used in second quantization the sign changes under $\Theta$ would suggest a transformation of $a_r \rightarrow a_r$ and $a_k \rightarrow a_{-k}$. If you are worried about the fact that $k$ represents a crystal momentum and not a true momentum you can just take the position transformation, which is perhaps more trustworthy, and use $a_k = \sum_r \, a_r \,\mathrm{exp}[ -\mathrm{i} k\cdot r]$ to verify $a_k \rightarrow a_{-k}$ directly.
Using these transformations you should be able to verify that the tight-binding Hamiltonian is invariant under time reversal in position and momentum space for a lattice with or without a basis. Keep in mind that you would generally take the complex conjugate of the coefficients in $H$, however in your case $t$ and $\epsilon_k$ are both real. Its important to remember though, mostly to make sure $H$ stays hermitian.
As far as your comment about $\sigma_y$, this is only necessary if you include spin. Spin changes sign under time reversal so $\Theta\,\mathbf{S}\,\Theta^\dagger = -\mathbf{S}$. In this case, we can formally write $\Theta = \mathrm{exp}[-i \pi J_y]\,K$, which is probably the relation you are alluding to.
According to J.J. Sakurai's Modern Quantum Mechanics one possible convention for the time-reversed angular momentum states is $\Theta | j,m\rangle = (-1)^m |j,-m\rangle$. This suggests that with spin indices the creation/annihilation operators transform like $a_{r,m} \rightarrow (-1)^m\,a_{r,-m}$ and $a_{k,m} \rightarrow (-1)^m \, a_{-k,-m}$ under time reversal. From what I understand, most spin Hamiltonians will be invariant under this transformation. An example when this is not the case would be in the presence of an external magnetic field which couples to the spins through a $\mathbf{S}\cdot \mathrm{B}-$like term.
It is interesting how even in the absence of an external field the groundstate of spin Hamiltonians can still sponanteously break the time reversal symmetry present in $H$, but rather than discuss this myself I will direct you to this very well written answer.
For each value of $k$, the above hamiltonian has energy contributions for an electron hopping from an $A$ to a $B$ site and vice versa. If we wanted to write the hamiltonian in matrix form (in a basis corresponding to $A$ and $B$ sites), it would look like
\begin{equation}
H_0 = \sum_{k,\alpha}\left(\begin{array}{c c}0& J_0\gamma_{k,\alpha}\\ J_0\gamma_{k,\alpha}^* & 0\end{array}\right)
\end{equation}
and we can diagonalize each term to get $\epsilon_{k,\alpha} = \pm J_0 |\gamma_{k,\alpha}|$.
Best Answer
You can try just working out the matrix elements between the different localised states, i.e. $a_1^{\dagger}\ldots a_m^{\dagger}|0\rangle$ etc. This may be rather tedious and require you to carefully keep track of phases due to the fermionic anti-symmetry. Probably a better way is to use the Jordan-Wigner transformation, which gives you an explicit matrix representation of the fermion operators. In 2D this representation of the Hamiltonian looks ugly, but it is still perfectly useful. You will be restricted to a fairly small lattice size, i.e. 10-20 sites depending on your computer and the efficiency of your code. It will be more efficient if you use the fermion number conservation to diagonalise within each number sector separately.
In the particularly simple case of a single particle on the lattice, the problem becomes rather easy, since exchange statistics are no longer relevant. In this case, you make the replacement $$a_{i,j}^{\dagger}a_{m,n} \to |i,j\rangle\langle m,n|,$$ where $$|i,j\rangle = a_{i,j}^{\dagger}|0\rangle.$$