Setup. The Laplacian $L = D - A$ has eigenvalues $0 = \mu_1 \leq \mu_2 \leq \cdots \leq \mu_n$ with eigenvector $\mathbf{1} = (1,\ldots,1)^T$ for $\mu_1 = 0$.

Step 1: Cofactor formula. We prove $t(G)$ equals any cofactor of $L$, say $\det(L_0)$ where $L_0$ is obtained by deleting the last row and column.

Cauchy-Binet formula. $\det(L_0) = \det(B_0 B_0^T)$ where $B$ is the oriented incidence matrix ($B_{ve} = 1$ if $v$ is the head of edge $e$, $-1$ if tail, 0 otherwise) and $B_0$ is $B$ with the last row deleted. By Cauchy-Binet: $\det(B_0 B_0^T) = \sum_{|S|=n-1} (\det B_0[S])^2$ where $S$ ranges over $(n-1)$-element subsets of edges.

Step 2: $\det B_0[S] \in \{-1, 0, 1\}$. The submatrix $B_0[S]$ corresponds to a subgraph $H_S$ with edge set $S$. If $H_S$ is a spanning tree: $B_0[S]$ is the reduced incidence matrix of a tree, which has $\det = \pm 1$ (proved by induction: a tree has a leaf, and expanding along the leaf's row gives a $\pm 1$ times a smaller tree incidence determinant). If $H_S$ is not a spanning tree: it contains a cycle or is disconnected, making $B_0[S]$ singular, so $\det = 0$.

Step 3: Therefore $\det(L_0) = \sum_T 1 = t(G)$, summing over spanning trees $T$.

Step 4: Spectral formula. Since $L$ has eigenvalues $0, \mu_2, \ldots, \mu_n$: $\det(L + tI) = t \prod_{i=2}^n (\mu_i + t)$ for all $t$. The cofactor $\det(L_0)$ equals $\frac{1}{n}\lim_{t\to 0}\frac{\det(L+tJ/n)}{1}$... More directly: by the matrix determinant lemma or direct calculation, any cofactor of $L$ equals $\frac{1}{n}\prod_{i=2}^n \mu_i$.

To see this: $L = Q\Lambda Q^T$ with $Q$ orthogonal and $\Lambda = \mathrm{diag}(0, \mu_2, \ldots, \mu_n)$. The cofactor matrix is $\mathrm{adj}(L) = \det(L)(L^{-1})$ if $L$ were invertible. Since $\mathrm{rank}(L) = n-1$, $\mathrm{adj}(L) = \prod_{i=2}^n \mu_i \cdot \frac{\mathbf{1}\mathbf{1}^T}{n}$ (the projection onto the null space, scaled). Every entry of $\mathrm{adj}(L)$ equals $\frac{1}{n}\prod_{i=2}^n\mu_i$, and each entry is a cofactor. $\square$