I guess mike's comment already answered your questions, but just in case, here's some additional information based on a more general framework. I'll consider your questions in the sequence: (2), (1), (3).
First, an observation: For continuous local martingales $M$ and $N$ with initial value zero, the process $\langle M, N\rangle$ is the unique (up to indistinguishability) continuous process with initial value zero and paths of finite variation such that $M_tN_t - \langle M,N\rangle_t$ is a local martingale. For the details of this, see for example Lemma IV.30.6 of Rogers & Williams: Diffusions, Markov Processes and Martingales, Volume 2.
(2). That $\int_0^t H_s d\langle M,N\rangle_s$ is continuous and of finite variation is essentially an analysis result and not a probabilistic result. Consider a continuous deterministic function $f:[0,\infty)\to\mathbb{R}$ with initial value zero and finite variation. This function $f$ induces for each $t\ge0$ a unique signed measure $\mu_t$ characterised by that $\mu_t((a,b]) = f(b)-f(a)$. Given for example a bounded Borel measurable function $g:[0,\infty)\to\mathbb{R}$, the integral $\int_0^t g(s)df(s)$ is defined as the Lebesgue integral of $g$ over $[0,t]$ with respect to $\mu_t$. Now, the signed measure $\mu_t$ has a Jordan-Hahn decomposition $\mu_t = \mu^+_t - \mu^-_t$ where $\mu^+$ and $\mu^-$ are nonnegative measures. We then obtain
$\int_0^t g(s)df(s) = \int_0^t g(s)d\mu^+_t(s) - \int_0^t g(s)d\mu^-_t(s)$.
In the case where $g$ is nonnegative, this shows that $\int_0^t g(s)df(s)$ as a function of $t$ is the difference between two increasing functions. Therefore, it is of finite variation. For the general case where $g$ is not nonnegative, write $g(s) = \max\{0,g(s)\} - \max\{0,-g(s)\}$, and the result follows for this case as well. We conclude that $\int_0^t g(s)df(s)$ has finite variation.
Furthermore, for any $t$, we have
$\int_0^t g(s)df(s) - \int_0^{t-}g(s)df(s) = \int 1_{\{t\}}g d\mu_t =g(t)\mu_t\{t\}=g(t)(f(t)-f(t-))=0$,
so $\int_0^t g(s)df(s)$ is continuous as a function of $t$ as well. Arguably, details about the definition and properties of integrals with respect to mappings of finite variation are left out of many expositions. For a few details on this, see Section IV.18 of Rogers & Williams: Diffusions, Markov processes and Martingales, Volume 2.
(1). I assume that $L_t = \int_0^t H_sdM_s$. By (2), the process $\int_0^t H_sd\langle M,N\rangle_s$ is continuous and has paths of finite variation. It also has initial value zero. Therefore, by the observation made earlier, in order to show $\langle L,N\rangle = \int_0^t H_sd\langle M,N\rangle_s$, it suffices to show that $L_tN_t - \int_0^t H_sd\langle M,N\rangle_s$ is a local martingale. The sign of $H$ is not relevant, this is taken care of in the uniqueness characterization of the quadratic covariation.
(3). First note that applying (1) with $N_t = \int_0^t H_sdM_s$, we obtain
$\langle \int_0^\cdot H_sdM_s\rangle_t =\langle \int_0^\cdot H_sdM_s,\int_0^\cdot H_sdM_s\rangle_t =\int_0^t H_sd\langle M,\int_0^\cdot H_udM_u\rangle_s$.
Next, applying (1) again, we obtain
$\langle M,\int_0^\cdot H_udM_u\rangle_s = \langle \int_0^\cdot H_udM_u,M\rangle_s = \int_0^s H_ud\langle M,M\rangle_u = \int_0^s H_ud\langle M\rangle_u$,
so integration with respect to $\langle M,\int_0^\cdot H_udM_u\rangle_s$ corresponds to integration with respect to the signed measure with density $H$ with respect to the measure induced by $\langle M\rangle$. Therefore,
$\int_0^t H_sd\langle M,\int_0^\cdot H_udM_u\rangle_s =\int_0^t H_s^2d\langle M\rangle_s$,
by the properties of the Lebesgue integral with respect to a measure with density. Collecting our conclusions, we obtain $\int_0^t H_s^2d\langle M\rangle_s = \langle \int_0^\cdot H_s dM_s\rangle_t$.