Does the Gauss-Bonnet theorem apply to non-orientable surfaces?
Orientability is not needed. Indeed, one can deduce unorientable version of Gauss-Bonnet theorem from the orientable one :
Given a (compact) nonorientable surface, say $M$, with metric $g$, consider its orientable double cover $\widetilde{M}$. The metric $g$ is naturally pulled back to a metric $\widetilde{g}$ on $\widetilde{M}$; i.e. locally, $\widetilde{M}$ is isometric to $M$. Then it is easy to see that the usual Gauss-Bonnet theorem on the orientable double cover implies the Gauss-Bonnet theorem on the nonorientable surface, because $2\pi\chi(M) = \pi\chi(\widetilde{M}) = \frac{1}{2}\int_{\widetilde{M}}{K} = \int_{M}K$, where $K$ denotes the Gaussian curvature.
Added :
One might want a version of Gauss-Bonnet theorem for surfaces with boundary. Indeed, the argument above can be applied to nonorientable surfaces with boundary. Let $M$ be a nonorientable surface with boundary $\partial M$. Then correspondingly its has an orientable double cover $\widetilde{M}$ with boundary $\partial \widetilde{M}$. As before, this double cover is locally isometric to $M$. Following the notation of Wikipedia, we have the Gauss-Bonnet theorem on $\widetilde{M}$ : $$\int_{\widetilde{M}}K + \int_{\partial\widetilde{M}}k_g = 2\pi \chi(\widetilde{M})$$ Now observe that each term are twice the corresponding term for $M$. In particular, $\int_{\partial\widetilde{M}}k_g = 2\int_{\partial M}k_g$ just because it is a 2-to-1 locally isometric double covering. As a result, we get $$\int_{M}K + \int_{\partial M}k_g = 2\pi \chi(M)$$ There is nothing special for nonorientable surfaces!
Added again :
Let me give you an example, a Möbius strip. Its Euler Characteristic is $0$. The most convenient metric on the Möbius strip is a flat metric; such Möbius strip can be realized as a quotient of a flat strip with parallel geodesic boundaries. With such a metric, $K=0$ and $k_g=0$, so the LHS of the Gauss-Bonnet theorem is $0$ as expected.
In general, for any metric on the Möbius strip, the Gauss-Bonnet theorem should hold just as well, because the LHS remains constant under smooth deformation of metric. (Any deformation is a composition of local deformations, and for local deformations, it is a corollary of the orientable version of Gauss-Bonnet theorem.)
Another easy proof is to cut your surface into small orientable pieces. It is not important whether your original surface is orientable or not. Gauss-Bonnet theorem holds for individual pieces, and when you glue them back, boundary terms corresponding to seams cancel out. In this way it is easy to see that orientability is not important at all in Gauss-Bonnet theorem for surfaces.
Following the OP's request, I'm posting the details of my computations for the explicit embedding of the Möbius strip $M$ in $\Bbb R^3$. The orientation of (halves) of the boundary curve turns out to be the crucial matter, as I'd suspected. Although the theoretical arguments are compelling, I remain confused about what's wrong with the following computations. I can now vouch that the numerics are correct.
Consider the parametrization $$x(u,v) = \big((2+v\sin(u/2))\cos u,(2+v\sin(u/2))\sin u,v\cos(u/2)\big), \quad 0\le u\le 2\pi, -1\le v\le 1.$$ Note that $B = x(0,1) = x(2\pi,-1)$ and $C=x(0,-1)=x(0,1)$. As you can check, this is an orthogonal parametrization, and the first fundamental form has coefficients $E = \|x_u\|^2 = 4+\frac34 v^2 - \frac12 v^2\cos u + 4v\sin(u/2)$ and $G=\|x_v\|^2 = 1$. I'm now going to use the standard formulas for $K$ and $\kappa_g$ in an orthogonal parametrization (see, e.g., my text, pp. 60 and 81): \begin{align*} K &= -\frac1{2\sqrt{EG}}\left(\Big(\frac{E_v}{\sqrt{EG}}\Big)_v + \Big(\frac{G_u}{\sqrt{EG}}\Big)_u\right) \\ \kappa_g &= \frac1{2\sqrt{EG}}\big({-}E_v u'(s) + G_u v'(s)\big)+\theta'(s), \end{align*} the latter for an arclength parametrization of the curve. Here $\theta$ is the angle the curve makes with $x_u$ at each point. In our case $\theta'=0$ everywhere. Note, also, that (forgetting about orientation issues for the moment) \begin{align*} \int_M K\,dA &= \iint_{[0,2\pi]\times [-1,1]} -\frac1{2\sqrt{EG}}\left(\Big(\frac{E_v}{\sqrt{EG}}\Big)_v + \Big(\frac{G_u}{\sqrt{EG}}\Big)_u\right)\underbrace{\sqrt{EG}\,du\,dv}_{dA} \\ &= -\frac12\int_0^{2\pi}\int_{-1}^1 \Big(\frac{E_v}{\sqrt E}\Big)_v\,dv\,du \\ &= -\frac12\int_0^{2\pi} \Big(\frac{E_v}{\sqrt E}\Big)\Big|_{v=1} - \Big(\frac{E_v}{\sqrt E}\Big)\Big|_{v=-1} \,du. \end{align*} Regarding the orientation issue, surely we'll agree that $|dA|$ should agree with this $dA$ if we remove one ruling of the surface (say the ruling from $B$ to $C$). So this integral should be the density integral.
Now, $E_v/\sqrt E$ is rather a mess, but, using Mathematica to do the numerical integration, we find that this integral is (approximately) $-1.97$. [The most basic check is that it's negative, as we have a non-developable ruled surface.]
We can now use Mathematica to evaluate the geodesic curvature integrals. We note that $ds = \frac{ds}{du}du$, so $\kappa_g\,ds = \big({-}\frac12 E_v/\sqrt E\big)u'(s)\,ds =\big({-}\frac12 E_v/\sqrt E\big)\,du$. On half the boundary circle, going from $B$ to $C$, $$-\frac12\int_0^{2\pi} \frac{E_v}{\sqrt E}\Big|_{v=1}du \approx -4.53,$$ and on the other half, going from $C$ to $B$, we have $$-\frac12\int_0^{2\pi} \frac{E_v}{\sqrt E}\Big|_{v=-1}du \approx -2.56.$$ So $\int_M K\,dA + \int_{\partial M} \kappa_g\,ds \approx -9.06$ (certainly not $0$, nor, indeed, an integer multiple of $2\pi$).
To double-check this computation, let's remove a tiny bit of our Möbius strip, say the region corresponding to $0\le u\le \varepsilon$. This leaves us an oriented surface, for sure. Its boundary has two extra pieces, $u=0$ (oriented downward) and $u=\varepsilon$ (oriented upward); these have no contribution, regardless of $\varepsilon$, since the $v$-curves are line segments and have no geodesic curvature. The main discrepancy, however, is the reversed orientation on the segment $v=1$. Indeed, for $\varepsilon$ very small, the $v=1$ integral is now approximately $+4.53$, and — mirabile dictu — note that $$-1.97 + 4.53 - 2.56 = 0,$$ as it should! (The exterior angles of the "rectangle" contribute the $2\pi\chi = 2\pi$.) But I emphasize that when we've removed a bit of the Möbius strip to make an oriented creature, we do not have (almost) the same boundary curve as the Möbius strip. This difference, as far as I'm concerned, is what messes up the Gauss-Bonnet Theorem. One more comment: The definition of $\kappa_g$ (as $\kappa\mathbf N\cdot(\mathbf n\times\mathbf T)$, where $\mathbf T,\mathbf N$ are the Frenet frame of the curve and $\mathbf n$ is the surface normal) makes it clear that when we have an oriented surface with boundary, the sign of $\kappa_g$ does not change if we reverse the orientation of the surface; for when we do, we change both $\mathbf n$ and $\mathbf T$ by a sign. But, if we interpreted $ds$ in this integral as a measure, as @SunghyukPark suggests, it surely must be consistent along the boundary circle of the Möbius strip, so we cannot just switch the sign of half the line integral. ...
EDIT: OK, I believe I've figured it out, much to my chagrin. We need to think about geodesic curvature intrinsically (as one does in a fancier proof of the Gauss-Bonnet Theorem). If $e_1$ is the unit tangent vector along $\partial M$ and $e_2$ is the inward-pointing normal to $\partial M$ in $M$, then, by definition, $\kappa_g = \nabla_{e_1}e_1\cdot e_2$. So, in fact, on the upper edge $v=1$ of our parametrizing rectangle, the formula for $\kappa_g$ we used earlier is off by a sign, and correcting this is equivalent to reversing the orientation on that upper integral. The correct values are $$\int_M K|dA| = -1.97 \qquad\text{and}\qquad \int_{\partial M}\kappa_g ds = +4.53-2.56 = +1.97.$$ So the sum is, in fact, $2\pi\chi(M) = 0$, as desired. I feel better now. :)
Here's the graph of the geodesic curvature (upper upper, lower lower):
This is more of an extended comment than a complete answer, but hope that it will close the question.
Applying a little Googlomagic (namely, searching for "gauss-bonnet non-orientable") it is possible to find out the following:
a paper by R. Palais's A Topological Gauss-Bonnet Theorem, J.Diff.Geom. 13 (1978) 385-398, where he mentions in passing that the Gauss-Bonnet theorem is easily generalized to the non-orientable case by considering measures.
an answer to this question with a feasible proof of the Gauss-Bonnet for the non-orientable case;
and many more interesting things, of course :)
On a side note, the Pfaffian has nothing to do with the orientability, bit rather with the dimension: it is defined in even dimensions (thus, in dimension 2 as well). See, maybe, here for the details.