I think I have figured it out.
We are looking for the topological response, hence it is enough to look at the response of any system described by a Hamiltonian adiabatically connected to the one we started with. In particular, we can use one with a flat spectrum: $H=\int\frac{d^2p}{(2\pi)^2}\psi^{*}(p) H(p)\psi(p)$ in which the relation $H(p)=U(p)\text{diag}(-I_{k},I_{n-k})U^{*}(p)$ holds locally. Here $k$ is the number of occupied bands and $U(k)=[v_1(p),...,v_{n}(p)]$ is an orthogonal matrix of eigenvectors of $H(p)$. In particular, the $k\times n$ matrix $S(p)=[v_1(p),...,v_{k}(p)]$ describes, locally, the occupied bundle. We can also write $H(p)=P^{\perp}(p)-P(p)$, where $P(p)=S(p)S^{*}(p)$ is the projector onto the fiber of the occupied bundle over $p$. The curvature of the bundle can be described as an endomorphism of the trivial bundle $\mathbb{T}^2\times \mathbb{C}^n$ by the expression $\Omega=PdP\wedge dPP=dP\wedge P^{\perp}dP$. The later expression will be useful in computing the effective action.
The idea is that we start with $Z_0=\text{Det}(G_0^{-1})$ and then perform the minimal coupling to an external gauge field. In a phase space representation ($(p,x)$-representation), and assuming that the external gauge field varies in scales larger than the system's typical scale, we can write the inverse of the new Green's function as $G^{-1}(p,x)\approx G_0^{-1}(p) -e A_{\mu}(x)\partial G_0^{-1}/\partial p_ \mu(p)$. Writing $\Sigma(p,x)=-eA_{\mu}(x)\partial G_0^{-1}/\partial p_ \mu(p)$ we have $\text{Det}(G^{-1})\approx Z_0\text{Det}(I+G_0\Sigma)$. The effective action is then obtained by the formal expansion $\log(\text{Det}(I+G_0\Sigma))=\text{Tr}\log(I+G_0\Sigma)\approx \text{Tr}(G_0\Sigma)-\frac{1}{2}\text{Tr}(G_0\Sigma G_0\Sigma)$. The partition function should be Gauge invariant, hence, we will only look for such terms. In $2+1$ dimensions, besides the Maxwell term, we are allowed to have the Chern-Simons term $S_{CS}(A)=(1/4\pi)\int A\wedge dA$. We are looking for the quadratic terms (hence, in the second term of the expansion written before) which have $A_\mu(x)$ and $\partial_\mu A_\nu(x)$. The products under the functional trace are in fact convolutions and, when we perform the transformation to the mixed position-momentum representation, we get a twisted product expansion, namely the Moyal product: $\int d^3x_2 A(x_1,x_2)B(x_2,x_3)\rightarrow A(p,x)B(p,x)+(i/2)\{A,B\}_{\text{PB}}(p,x)+...$, where $\{.,.\}_{\text{PB}}$ is the Poisson Bracket. The first two terms of the expansion are enough for us. If we look for the mentioned terms we end up with the contribution
$(ie^{2}/2)(1/3!)(\int d^3p/(2\pi)^3 \text{tr} (G_0\partial G_0^{-1}/\partial p_\mu G_0\partial G_0^{-1}/\partial p_\nu G_0\partial G_0^{-1}/\partial p_\lambda) \varepsilon_{\mu\nu\lambda}) \int A\wedge dA$. By performing the integral over the frequency $p_0$, the later can be shown to be equal to $2\pi i\sigma_{H} S_{CS}(A)$, with $\sigma_{H}=(e^{2}/2\pi)\times \int_{\mathbb{T}^2}\text{tr}(i\Omega/2\pi)\equiv (e^{2}/2\pi)\times c_1$ ($c_1$ denotes the first Chern number of the occupied bundle), as expected.