How can off-diagonal Hamiltonian matrix elements be nonzero when the overlap matrix is diagonal?
I believe the confusion lies in some ambiguity in the use of the term overlap.

For the $S_{AB}=0$ approximation of Hückel theory, overlap is used in the specific mathematical sense of the "overlap integral": 

$$
S_{AB} = \left&lt;\varphi_A|\varphi_B\right&gt; = 
\int_{\vec r}{\varphi_A^*\!\left(\vec r\right)\varphi_B\!\left(\vec r\right)\mathrm{d}\vec r}
$$

However, this does not mean that at each point in $\mathbb{R}^3$ only one of $\varphi_A$ or $\varphi_B$ is allowed to bear a non-zero value.  In fact, at most points in $\mathbb R^3$, $\varphi_A$ and $\varphi_B$ will both bear non-zero values. Since $\varphi_A$ and $\varphi_B$ will in general have "co-located regions of non-zero values," the two functions can be described qualitatively as overlapping, in the same sense as in your second quote. It is the overlap of the functions in this more general, qualitative sense, that allows $H_{AB}=\left&lt;\varphi_A^*\middle|\hat H\middle|\varphi_B\right&gt;\neq0$.

(To note,  in order for Hückel's approximation to be exact, the orbitals would have to be selected such that the specific integral $S_{AB}$ over all $\mathbb R^3$ equals zero$^\ddagger$. I believe Hückel did not actually do this, but instead just coerced to zero any nonzero $S_{AB}$.)

To help explain more explicitly why this is, I'll first write out the integral implied by the Dirac notation for $H_{AB}$:

$$
H_{AB}=\left&lt;\varphi_A^*\middle|\hat H\middle|\varphi_B\right&gt;=
\int_{\vec r}{\varphi_A^*\!\left(\vec r\right)\hat H\!\left[\varphi_B\!\left(\vec r\right)\right]\mathrm{d}\vec r}
$$

Then, I'll introduce two (real-valued) one-dimensional functions defined on the closed domain $x=[0,1]$ to use for demonstration purposes:

$$
\varphi_A\!\left(x\right) = \sin\!\left(2\pi x\right) \\
\varphi_B\!\left(x\right) = \sin\!\left(3\pi x\right)
$$

$\hspace{10mm}$

It is straightforward to show that $S_{AB}=0$ for these two functions.  One way is to plug the overlap integral into Wolfram Alpha. Another option is to do it by hand.  In the below figure, the green curve is a plot of the product $\varphi_A\!\left(x\right)\cdot\varphi_B\!\left(x\right)$, and the orange curve is a plot of the cumulative integral of this product:

$$
\int_0^x{\varphi_A\!\left(x'\right)\cdot\varphi_B\!\left(x'\right)\mathrm{d}x'}
$$

$\hspace{10mm}$

As can be seen, the integral over the entire domain is zero, satisfying $S_{AB}=0$.

Now, let's introduce a simple "Hamiltonian"$^\dagger$:

$$
\hat H = x
$$

The next figure is a replicate of the first one above, with $\hat H\!\left[\varphi_B\right]$ also plotted as the dashed red curve:

$\hspace{10mm}$

As in the second figure above, one last plot shows the local value of $\varphi_A\!\left(x\right)\hat H\!\left[\varphi_B\!\left(x\right)\right]$ along with the cumulative integral:

$$
\int_0^x{\varphi_A\!\left(x'\right)\cdot\left[x'\cdot\varphi_B\!\left(x'\right)\right]\mathrm{d}x'}
$$

$\hspace{10mm}$

The value of $H_{AB}$ here is clearly not equal to zero.  (Wolfram Alpha concurs.) Thus, it is absolutely possible to have $H_{AB} \neq S_{AB} = 0$.  

To note, the principle exhibited here on a closed domain in $\mathbb R^1$ applies equally well to orbitals on an infinite domain in $\mathbb R^3$. Note also, if I were actually performing a quantum mechanical calculation here, I would have been a lot more careful about normalizing my wavefunctions (viz., ensuring $S_{ii}=1$).

$^\ddagger$Numerous methods are known for generating orthogonalized orbitals with $S_{AB}=0$ strictly, from a general set of orbitals where $S_{AB}\neq0$ occurs for some non-null subset of the $\left(\!A,\!B\right)$. To my knowledge, the most commonly used orthogonalization for atomic orbitals is Löwdin's method (see this PDF). The Gram-Schmidt orthogonalization (usually in its "modified" form) is used in a variety of other areas of quantum chemistry. Wikipedia also has a pretty solid article on orthogonalization methods in general.

$^\dagger$You'd never use this as a Hamiltonian operator. Further, I'm not sure if it technically even is a Hamiltonian operator. (I think it'd work as a dipole moment operator?) But, whatever, leave me alone, it's for purposes of illustration.