It looks like magic because the proof of Kirchhoff's Matrix Tree Theorem is nontrivial.
It relies on several algebraic properties of the matrix constructed in steps 1-3, which is called the Laplacian matrix of the graph.
Let $A$ be the adjacency matrix, and let $D$ be the diagonal matrix with the degrees of the nodes on the diagonal. Steps 1-3 build the matrix $L = D-A$: this is the Laplacian matrix of the graph.
So the question becomes, why should a minor of the Laplacian matrix be equal to the number of spanning trees of the graph?
To answer this question, we need to use other properties of $L$. One property is that we can also write $L$ as $B \cdot B^\top$ where $B$ is the oriented incidence matrix of the graph, that is, a matrix with one row per node $u$ and one column per edge $e$, where entry $B_{ue}=1$ if $e=\{u,v\}$ with $u<v$, $B_{ue}=-1$ if $e=\{u,v\}$ with $v < u$, and $B_{ue}=0$ if $e$ does not contain $u$.
For example, for a triangle graph, the Laplacian matrix is
$$ L=\left[\begin{array}{ccc} 2 & -1 & -1 \\ -1 & 2 & -1 \\ -1 & -1 & 2 \end{array}\right] $$
and the incidence matrix is
$$ B=\left[\begin{array}{ccc} 1 & 1 & 0 \\ -1 & 0 & 1 \\ 0 & -1 & -1 \end{array}\right]. $$
Note that the inner product of two distinct rows of $B$ is $-1$ if the corresponding nodes are adjacent, and 0 if they are not. The inner product of a row of $B$ with itself equals the degree of the corresponding node. (Just use the definition of $B$.) This is why $L=B\cdot B^\top$.
Assuming the graph is connected, the rank of $B$ is exactly $n-1$ where $n$ is the number of nodes of the graph. It is not $n$ because if you sum all rows, you get a zero vector. It is $n-1$ because the number of independent columns of $B$ is $n-1$: if you take any subset of columns corresponding to the set of edges of any spanning tree of the graph, these columns will be independent; the only way to have a dependent set of columns is when the edges selected contain a cycle.
If we delete a row and corresponding column of $L$, we obtain a matrix $M$. If we delete the corresponding row of $B$, we obtain a matrix $C$. Now $M=C \cdot C^\top$ for the same reasons that $L=B \cdot B^\top$. We are interested in $\det(M)=\det(CC^\top)$. In the example of the triangle graph above,
$$ M=\left[\begin{array}{cc} 2 & -1 \\ -1 & 2 \end{array}\right] $$
and
$$ C=\left[\begin{array}{ccc} 1 & 1 & 0 \\ -1 & 0 & 1 \end{array}\right]. $$
The last step is to invoke another algebraic fact, the Cauchy-Binet formula. This says that for any matrix $C$,
$$ \det(C C^\top) = \det (C_{K_1}) \det(C_{K_1}^\top) + \ldots + \det (C_{K_q}) \det(C_{K_q}^\top) $$ where each $C_{K_i}$ denotes an $(n-1) \times (n-1)$ submatrix of $C$ and the sum runs over all such submatrices. Thus,
$$ \det(C C^\top) = \det (C_{K_1})^2 + \ldots + \det (C_{K_q})^2. $$
However, $\det(C_{K_i})=\pm 1$ if the edges corresponding to $C_{K_i}$ comprise a spanning tree, while $\det(C_{K_i})=0$ otherwise (because then the set contains a cycle and the corresponding columns are not independent). Thus, the value of the sum is exactly the number of spanning trees.
