Here is how the heat kernel proof of Atiyah-Singer goes at a high level. Let $(\partial_t - \Delta)u = 0$ and define the heat kernel (HK) or Green function via $\exp(-t\Delta):u(0,\cdot) \rightarrow u(t,\cdot)$. The HK derives from the solution of the heat equation on the circle:
$u(t,\theta) = \sum_n a_n(t) \exp(in\theta) \implies a_n(t) = a_n(0)\cdot \exp(-tn^2)$
For a sufficiently nice case the solution of the heat equation is $u(t,\cdot) = \exp(-t\Delta) * u(0,\cdot)$.
The hard part is building the HK: we have to compute the eigenstuff of $\Delta$ (this is the Hodge theorem). But once we do that, a miracle occurs and we get the
Atiyah-Singer Theorem: The supertrace
of the HK on forms is constant: viz.
$\mathrm{Tr}_s \exp(-t\Delta) = \sum_k (-1)^k \,\mathrm{Tr} \exp(-t\Delta^k) = \mathrm{const}.$
For $t$ large, this can be evaluated topologically; for small $t$, it can be evaluated analytically as an integral of a characteristic class.
Edit per Qiaochu's clarification
This article of Kotake (really in here as the books seem to be mixed up) proves Riemann-Roch directly using the heat kernel.
When $E\to M$ is an oriented vector bundle of rank $2n$ over a
compact manifold $M$, it has a well-defined de Rham Euler class $e(E)$
in $H^{2n}_{dR}(M)$, and a representative $2n$-form for $e(E)$
can be computed as follows:
Fix a positive definite inner product $\langle,\rangle$ on $E$. (Since any two such inner products are equivalent under automorphisms of $E$, it doesn't matter which one.) Let $\nabla$ be a $\langle,\rangle$-orthogonal connection on $E$, and let $K^\nabla$ denote the curvature of $\nabla$, regarded as a $2$-form with values in ${\frak{so}}(E)$. Let $e(\nabla) = \mathrm{Pf}(K^\nabla/2\pi)$, which is a well-defined $2n$-form on $M$. (N.B.: The definition of $\mathrm{Pf}$ requires both the inner product and the orientation of $E$.) Then $e(E)=\bigl[e(\nabla)\bigr]\in H^{2n}_{dR}(M)$. (That $\bigl[e(\nabla)\bigr]$ is independent of the choice of $\nabla$
is one of the first results proved in Chern-Weil theory.)
If $M$ is a compact, oriented $2n$-manifold, then the value of $e(E)$ on $[M]$, the fundamental class of $M$, can be computed as follows: Let $Y$ be a section of $E$ that has only a finite number of zeroes. (By Whitney transversality, the generic section of $E$ satisfies this condition.) Using the orientations of $E$ and $M$, one defines the index
of an isolated zero $z$ of $Y$, which is an integer $\iota_Y(z)$. Then
$$
e(E)\bigl([M]\bigr) = \int_M e(\nabla) = \sum_{z\in Z}\ \iota_Y(z).
$$
By the Poincaré-Hopf Theorem, when $E = TM$ (as oriented bundles),
this sum is equal to the Euler characteristic of $M$, which explains
why $e(E)$ is called the Euler class of $E$.
To prove this result, one chooses an orthogonal connection $\nabla$ on $E$ that depends on $Y$ and whose Euler form $e(\nabla)$ is easy to evaluate explicitly. This can be done as follows:
Let $Z\subset M$ be the (finite) zero set of $Y$ and let $U\subset M$ be an open neighborhood of $Z$ that consists of disjoint, smoothly embedded $2n$-balls, one around each element of $Z$. Let $\phi$ be a smooth function on $M$ that is identically equal to $1$ on an open neighborhood of $Z$ and whose support $K\subset U$ is a disjoint union of closed, smoothly embedded balls,
one around each element of $Z$.
Now, $E$ is trivial over $U$,
so choose a positively oriented $\langle,\rangle$-orthonormal basis of sections
$s_1,\ldots, s_{2n}$ of $E$ over $U$ and let $\nabla_1$ be the (flat) connection on $E_U$ for which these sections are parallel. Let $\bar Y = Y/|Y|$ be the normalized unit section defined on $M\setminus Z$, and define a second connection on $U$ by
$$
\nabla_2 s = \nabla_1 s
+ (1{-}\phi)\bigl( \bar Y\otimes \langle s,\nabla_1 \bar Y\ \rangle
- \langle s,\bar Y\ \rangle\ \nabla_1\bar Y\ \bigr).
$$
It is easily verified that this formula does define a connection on $U$,
that $\nabla_2$ is $\langle,\rangle$-orthogonal, and that,
outside of $K$ (i.e., where $\phi\equiv0$), the vector field $\bar Y$ is
$\nabla_2$-parallel. (Note that, because $\phi\equiv1$ near $Z$, where $Y$ is
not defined, this formula does extend smoothly across $Z$, agreeing with $\nabla_1$
on a neighborhood of $Z$.)
Meanwhile, on $M\setminus K$, write $E$ as a $\langle,\rangle$-orthogonal direct sum $E = \mathbb{R}\cdot Y \oplus E'$ and choose a $\langle,\rangle$-compatible connection $\nabla_3$ on $E$ over $M\setminus K$ that preserves this splitting. Using a partition of unity subordinate to the open cover of $M$ defined by $U$ and $M\setminus K$, construct a $\langle,\rangle$-orthogonal connection $\nabla$ that agrees with $\nabla_2$ on $K$ and with $\nabla_3$ on $M\setminus U$, and for which the line bundle $\mathbb{R}\cdot Y\subset E$
is parallel on $M\setminus K$.
Now, on $M\setminus K$, the curvature of $\nabla$ takes values in ${\frak{so}}(E')\subset{\frak{so}}(E)$, so $e(\nabla) = \textrm{Pf}\bigl(K^{\nabla}/(2\pi)\bigr)$ vanishes identically outside of $K$. It suffices now to evaluate the
integral of $e(\nabla)$ over a single component $B$ of $K$, which may be
assumed to be the unit ball in $\mathbb{R}^{2n}$, so restrict attention to
this case. Write $\bar Y = s_1 u_1 +\cdots + s_{2n}\ u_{2n}$, and note that,
on $B$, one has, by definition,
$$
\nabla s_i = \nabla_2 s_i
= \sum_{j=1}^{2n} s_j\otimes (1{-}\phi)(u_j\ du_i - u_i\ du_j),
$$
i.e., the connection $1$-forms of $\nabla$ in this basis are
$\omega_{ij} = (1{-}\phi)(u_i\ du_j - u_j\ du_i)$.
Using the identity ${u_1}^2 +\cdots + {u_{2n}}^2 = 1$,
the curvature forms are easily computed to satisfy
$$
\Omega_{ij} = d\omega_{ij} + \omega_{ik}\wedge\omega_{kj}
= (1{-}\phi^2)\ du_i\wedge du_j - d\phi\wedge(u_i\ du_j - u_j\ du_i).
$$
At this point, you have to know the definition of the Pfaffian. (I'll wait while you look it up.) Using the fact that the result has to be invariant under $SO(2n)$-rotations, it is easy to show that, on $B$, one has
$$
e(\nabla) = \textrm{Pf}\left(\frac{\Omega}{2\pi}\right)
= c_n (1{-}\phi^2)^{n-1} d\phi\wedge u^*(\Upsilon)
$$
for some universal constant $c_n$ and where $\Upsilon$ is the $SO(2n)$-invariant $2n$-form
of unit volume on $S^{2n-1}\subset\mathbb{R}^{2n}$ and $u: B\setminus z\to S^{2n-1}$
is $u = (u_1,\ldots, u_{2n})$. (I'll let you evaluate the constant $c_n$. This can be done in a number of ways, but it's essentially a combinatorial exercise.) In particular,
$$
e(\nabla) = d\bigl(P_n(1{-}\phi)\ u^*(\Upsilon)\bigr)
$$
where $P_n(t)$ is the polynomial in $t$ (of degree $2n{-}1)$) that vanishes at $t=0$ and
satisfies $P'_n(1{-}t) = -c_n(1{-}t^2)^{n-1}$. By Stokes' Theorem, one has
$$
\int_B e(\nabla) = P_n(1) \int_{\partial B}u^*(\Upsilon)
= P_n(1)\ \textrm{deg}(u:\partial B\to S^{2n-1}) = P_n(1)\ \iota_Y(z)
$$
Thus, it follows that there is a universal constant $C_n = P_n(1)$ such that
$$
e(E)\bigl([M]\bigr) = \int_M e(\nabla) = C_n \ \sum_{z\in Z}\ \iota_Y(z).
$$
Now, $C_1 = 1$, as one can show by computing with the constant curvature metric on the $2$-sphere and noting that, for $Y$ the gradient of the height function on $S^2$, the sum of the indices is $2$. Finally, by properties of the Pfaffian and the index, one sees that,
if such a formula as above is to hold, then one must have $C_{n+m}=C_nC_m$. Thus, $C_n=1$ for all $n$, and the formula is proved.
(Of course, one can avoid these final tricks for evaluating the constants by carefully doing the combinatorial exercise, but I'll leave that to the curious.)
Best Answer
The answer is yes, to both questions.
First question first. For any geodesic $n$-gon $P$ on $M$, i.e., a simply connected region of $M$ whose boundary consists of $n$-geodesic arcs, define
$$ \delta(P)= \mbox{sum of the angles of $T$}-(n-2)\pi. $$
Note that if $M$ were flat, then the defect $\delta(P)$ would be $0$. This quantity has a remarkable property: if $P= P'\cup P''$, where $P, P', P''$ are geodesic polygons, then
$$ \delta(P)=\delta(P')+\delta(P'')-\delta(P'\cap P'') $$
which shows that $\delta$ behaves like a finitely additive measure. It can be extended to a countably additive measure on $M$, and as such, it turns out to be absolutely continuous with respect to the the volume measure $dV_g$ defined by the Riemann metric $g$. $\newcommand{\bR}{\mathbb{R}}$ Thus we can find a function $\rho: M\to \bR$ such that
$$\delta= \rho dV_g. $$
More concretely for any $p\in M$ we have
$$\rho(p) =\lim_{P\searrow p} \frac{\delta(P)}{{\rm area}\;(P)}, $$
where the limit is taken over geodesic polygons $P$ that shrink down to the point $p$. In fact
$$ \rho(p) = K(p). $$
Now observe that if we have a geodesic triangulation$\newcommand{\eT}{\mathscr{T}}$ $\eT$ of $M$, the combinatorial Gauss-Bonnet formula reads
$$\sum_{T\in\eT} \delta(T)=2\pi \chi(M). $$
On the other hand
$$\delta(T) =\int_T \rho(p) dV_g(p), $$
and we deduce
$$ \int_M \rho(p) dV_g(g)=\sum_{T\in \eT}\int_T \rho(p) dV_g(p) =\sum_{T\in\eT} \delta(T)=2\pi \chi(M). $$
For more details see these notes for a talk I gave to first year grad students a while back.
As for the second question, perhaps the most general version of Gauss-Bonnet uses the concept of normal cycle introduced by Joseph Fu.
This is a rather tricky and technical subject, which has an intuitive description. Here is roughly the idea.
To each compact and reasonably behaved subset $S\subset \bR^n$ one can associate an $(n-1)$-dimensional current $\newcommand{\bN}{\boldsymbol{N}}$ $\bN^S$ that lives in $\Sigma T\bR^n =$ the unit sphere bundle of the tangent bundle of $\bR^n$. Think of $\bN^S$ as oriented $(n-1)$-dimensional submanifold of $\Sigma T\bR^n$. The term reasonably behaved is quite generous because it includes all of the examples that you can produce in finite time (Cantor-like sets are excluded). For example, any compact, semialgebraic set is reasonably behaved.
How does $\bN^S$ look? For example, if $S$ is a submanifold, then $\bN^S$ is the unit sphere bundle of the normal bundle of $S\hookrightarrow \bR^n$.
If $S$ is a compact domain of $\bR^n$ with $C^2$-boundary, then $\bN^S$, as a subset of $\bR^N\times S^{n-1}$ can be identified with the graph of the Gauss map of $\partial S$, i.e. the map
$$\bR^n\supset \partial S\ni p\mapsto \nu(p)\in S^{n-1}, $$
where $\nu(p)$ denotes the unit-outer-normal to $\partial S$ at $p$.
More generally, for any $ S$, consider the tube of radius $\newcommand{\ve}{{\varepsilon}}$ around $S$
$$S_\ve= \bigl\lbrace x\in\bR^n;\;\; {\rm dist}\;(x, S)\leq \ve\;\bigr\rbrace. $$
For $\ve $ sufficiently small, $S_\ve$ is a compact domain with $C^2$-boundary (here I'm winging it a bit) and we can define $\bN^{S_\ve}$ as before. $\bN^{S_\ve}$ converges in an appropriate way to $\bN^{S}$ as $\ve\to 0$ so that for $\ve$ small, $\bN^{S_\ve}$ is a good approximation for $\bN^S$. Intuitively, $\bN^S$ is the graph of a (possibly non existent) Gauss-map.
If $S$ is a convex polyhedron $\bN^S$ is easy to visualize. In general $\bN^S$ satisfies a remarkable additivity
$$\bN^{S_1\cup S_2}= \bN^{S_1}+\bN^{S_2}-\bN^{S_1\cap S_2}. $$
In particular this leads to quite detailed description for $\bN^S$ for a triangulated space $S$.
Where does the Gauss-Bonnet formula come in? As observed by J. Fu, there are some canonical, $O(n)$-invariant, degree $(n-1)$ differential forms on $\Sigma T\bR^n$, $\omega_0,\dotsc, \omega_{n-1}$ with lots of properties, one being that for any compact reasonable subset $S$
$$\chi(S)=\int_{\bN^S}\omega_0. $$
The last equality contains as special cases the two formulae you included in your question.
I am aware that the last explanations may feel opaque at a first go, so I suggests some easier, friendlier sources.
For the normal cycle of simplicial complexes try these notes. For an exposition of Bernig's elegant approach to normal cycles try these notes.
Even these "friendly" expositions with a minimal amount technicalities could be taxing since they assume familiarity with many concepts.
Last, but not least, you should have a look at these REU notes on these subject. While the normal cycle does not appear, its shadow is all over the place in these beautifully written notes.
Also, see these notes for a minicourse on this topic that I gave a while back.