We develop the twistor theory of $G$-structures for which
the (linear) Lie algebra of the structure group contains an
involution, instead of a complex structure. The twistor space
$Z$ of such a $G$-structure is endowed with a field of involutions
$\mathcal{J}\in \Gamma (\End TZ)$ and a $\mathcal{J}$-invariant
distribution $\mathcal{H}_{Z}$. We study the conditions for
the integrability of $\mathcal{J}$ and for the (para-)holomorphicity
of $\mathcal{H}_{Z}$. Then we apply this theory to para-quaternionic
Kähler manifolds of non-zero scalar curvature, which
admit two natural twistor spaces $(Z^{\epsilon},\mathcal{J},\mathcal{H}_{Z})$,
$\epsilon=\pm 1$, such that $\mathcal{J}^{2}=\epsilon \Id$.
We prove that in both cases $\mathcal{J}$ is integrable (recovering
results of Blair, Davidov and Mu\u{s}karov) and that $\mathcal{H}_{Z}$
defines a holomorphic ($\epsilon=-1$) or para-holomorphic
($\epsilon=+1$) contact structure. Furthermore, we determine
all the solutions of the Einstein equation for the canonical
one-parameter family of pseudo-Riemannian metrics on $Z^{\epsilon}$.
In particular, we find that there is a unique Kähler-Einstein
($\epsilon=-1$) or para-Kähler-Einstein ($\epsilon=+1$)
metric. Finally, we prove that any Kähler or para-Kähler
submanifold of a para-quaternionic Kähler manifold
is minimal and describe all such submanifolds in terms of
complex ($\epsilon=-1$), respectively, para-complex ($\epsilon=+1$)
submanifolds of $Z^{\epsilon}$ tangent to the contact distribution.