This is only an answer to the request for references about the expansion and an interpretation of $h_2$ from the standpoint of the metric $x^2g$.
Theorem 7.4 in the book The Ambient Metric proves that any hyperbolic metric on a conformally compact manifold has an expansion of the type you specify.
Note that in the three-dimensional case you are asking about, $h_4$ is determined by $h_2$, but $h_2$ isn’t determined solely by the conformal boundary.
Now let me explain how I think about $h_2$, starting in the general framework of Poincaré-Einstein manifolds.
Let $(M^{n+1},g_+)$, $n \geq 2$, be a Poincaré-Einstein manifold and let $x$ be a geodesic defining function.
Then
$$ g_+ = x^{-2}\bigl(dx^2 + h_0 + x^2h_2 + o(x^2) \bigr) $$
near $\partial M$.
Denote $g := x^2g_+$.
Let $P^{g_+}$ and $P^g$ denote the Schouten tensors of $g_+$ and $g$, respectively.
(Since $\dim M \geq 3$, the Schouten tensors are defined.)
By assumption, $P^{g_+} = -\frac{1}{2}g_+$.
Since $x$ is a geodesic defining function, the conformal transformation law for the Schouten tensor implies that
$$ P^g = -x^{-1}\nabla^2x , $$
where the Hessian is taken with respect to $g$.
Direct computation gives that
$$ \nabla^2x = -xh_2 + o(x) $$
near $\{ x = 0 \}$.
In particular, the restriction of $P^g$ to $\partial M$ is
$$ \tag{1}\label{1} P^g\rvert_{\partial M} = -h_2 . $$
Consider now the case of a hypersurface $\Sigma^n$, $n \geq 2$, in a Riemannian manifold $(M^{n+1},g)$.
Denote by $h := g\rvert_{\Sigma}$ the induced metric on $\Sigma$, by $L$ the second fundamental form, and by $H := \frac{1}{n}\operatorname{tr}_h L$ the mean curvature.
Set
$$ \tag{2}\label{2} \mathcal{P} := P\rvert_\Sigma + HL - \frac{1}{2}H^2h . $$
One can check that $\mathcal{P}$ transforms like a Schouten tensor; i.e.
$$ \mathcal{P}^{u^{-2}g} = \mathcal{P}^g + u^{-1}\nabla_h^2u - \frac{1}{2}\lvert\nabla_h u\rvert_h^2 h . $$
Note in particular that this depends only on $u\rvert_\Sigma$, and hence $\mathcal{P}$ depends only on the conformal class $[g]$ and a choice of metric $h \in [g\rvert_\Sigma]$.
The motivation for this definition comes from the Gauss equations: if $n \geq 3$, then
$$ \mathcal{P} = \overline{P} + F , $$
where $\overline{P}$ is the Schouten tensor of $h$ and $F$ is the Fialkow tensor (see Equation (3.34) here, for example).
Two key properties of $F$ are that it is conformally invariant and it vanishes on the boundary of any Poincar'e-Einstein manifold.
That is, $\mathcal{P}$ can be thought of as a generalization of the Schouten tensor of $h$ to the case $n=2$, where the latter is not defined.
In other words, $\mathcal{P}$ defines a Möbius structure on $\Sigma$ from the conformal structure on $(M^{n+1},g)$, and this coincides with the usual Möbius (=conformal) structure when $n \geq 3$ (this is discussed further in Section 3.5 here).
Finally, coming back to the case of Poincar'e-Einstein manifolds.
Equation \eqref{1} implies that $h_2 = -\mathcal{P}^g$, meaning that the tensor $\mathcal{P}$ above is exactly $P_{ij}$ appearing in Theorem 7.4 of The Ambient Metric.
Moreover, Equation \eqref{2} and the discussion surrounding $\mathcal{P}$ implies that one could compute $\mathcal{P}$ as the limit as $x \to 0^+$ of the same quantities computed with respect to $g_+$;
i.e. $\mathcal{P}$ can be computed by using the second fundamental form of the level sets $\{ x = \varepsilon \}$ and taking the limit $\varepsilon\to0^+$.
Finally, $\mathcal{P}$ determines the Möbius structure on $\partial M$, and hence which element of the Teichmüller space is being represented.
I don't know anything about McMullen's quasi-fuchsian reciprocity, so cannot help you there.
