[proofplan]
We prove the subgroup property directly from the defining behaviour of horizontal lifts. The constant loop gives the identity element. Reversing a loop reverses its horizontal transport, giving inverses after using equivariance of principal parallel transport under the right $G$-action. Concatenating loops corresponds to composing their parallel transports, and equivariance identifies the endpoint with right multiplication by the product in $G$.
[/proofplan]
[step:Fix the notation for principal parallel transport]
Let $x_0 := \pi(p)$. For every piecewise smooth loop
\begin{align*}
\gamma: [0,1] \to M
\end{align*}
based at $x_0$, define the principal parallel transport map along $\gamma$ by
\begin{align*}
\operatorname{PT}_{\gamma}: \pi^{-1}(x_0) &\to \pi^{-1}(x_0)
\end{align*}
where $\operatorname{PT}_{\gamma}(q)$ is the endpoint at time $1$ of the $\omega$-horizontal lift of $\gamma$ starting at $q \in \pi^{-1}(x_0)$.
By definition,
\begin{align*}
\operatorname{Hol}_p(\omega)=\{g \in G : \operatorname{PT}_{\gamma}(p)=p g \text{ for some piecewise smooth loop } \gamma \text{ based at } x_0\}.
\end{align*}
We use the standard equivariance property of principal parallel transport: for every loop $\gamma$ based at $x_0$, every $q \in \pi^{-1}(x_0)$, and every $a \in G$,
\begin{align*}
\operatorname{PT}_{\gamma}(q a)=\operatorname{PT}_{\gamma}(q)a.
\end{align*}
This holds because the connection is principal, so the horizontal distribution is invariant under the right $G$-action.
[/step]
[step:Use the constant loop to obtain the identity element]
Let $e \in G$ denote the identity element. Define the constant loop
\begin{align*}
c: [0,1] &\to M
\end{align*}
by $c(t)=x_0$ for all $t \in [0,1]$. The constant path
\begin{align*}
\widetilde{c}_p: [0,1] &\to P
\end{align*}
given by $\widetilde{c}_p(t)=p$ is horizontal and satisfies $\widetilde{c}_p(0)=p$. Hence $\operatorname{PT}_c(p)=p=p e$, so $e \in \operatorname{Hol}_p(\omega)$.
[/step]
[step:Reverse a loop to obtain the inverse holonomy element]
Let $g \in \operatorname{Hol}_p(\omega)$. Then there exists a piecewise smooth loop
\begin{align*}
\gamma: [0,1] &\to M
\end{align*}
based at $x_0$ such that $\operatorname{PT}_{\gamma}(p)=p g$.
Define the reversed loop
\begin{align*}
\gamma^{-1}: [0,1] &\to M
\end{align*}
by $\gamma^{-1}(t)=\gamma(1-t)$. If $\widetilde{\gamma}_p: [0,1] \to P$ is the horizontal lift of $\gamma$ starting at $p$, then the path
\begin{align*}
\eta: [0,1] &\to P
\end{align*}
defined by $\eta(t)=\widetilde{\gamma}_p(1-t)$ is the horizontal lift of $\gamma^{-1}$ starting at $p g$ and ending at $p$. Therefore
\begin{align*}
\operatorname{PT}_{\gamma^{-1}}(p g)=p.
\end{align*}
By equivariance of principal parallel transport,
\begin{align*}
\operatorname{PT}_{\gamma^{-1}}(p g)=\operatorname{PT}_{\gamma^{-1}}(p)g.
\end{align*}
Thus $\operatorname{PT}_{\gamma^{-1}}(p)g=p$, and since the right $G$-action on each fiber is free, this implies
\begin{align*}
\operatorname{PT}_{\gamma^{-1}}(p)=p g^{-1}.
\end{align*}
Hence $g^{-1} \in \operatorname{Hol}_p(\omega)$.
[guided]
Let $g \in \operatorname{Hol}_p(\omega)$. By the definition of holonomy, there is a piecewise smooth loop
\begin{align*}
\gamma: [0,1] &\to M
\end{align*}
based at $x_0=\pi(p)$ such that horizontal transport along $\gamma$ carries $p$ to $p g$:
\begin{align*}
\operatorname{PT}_{\gamma}(p)=p g.
\end{align*}
To produce the inverse element, we run the same base loop in the opposite direction. Define
\begin{align*}
\gamma^{-1}: [0,1] &\to M
\end{align*}
by $\gamma^{-1}(t)=\gamma(1-t)$. If $\widetilde{\gamma}_p: [0,1] \to P$ is the horizontal lift of $\gamma$ starting at $p$, then the path
\begin{align*}
\eta: [0,1] &\to P
\end{align*}
defined by $\eta(t)=\widetilde{\gamma}_p(1-t)$ projects to $\gamma^{-1}$, starts at $\widetilde{\gamma}_p(1)=p g$, and ends at $\widetilde{\gamma}_p(0)=p$. Reversing a horizontal path is still horizontal, because its tangent vectors remain in the same horizontal subspaces with their signs reversed. Hence $\eta$ is the horizontal lift of $\gamma^{-1}$ starting at $p g$, and so
\begin{align*}
\operatorname{PT}_{\gamma^{-1}}(p g)=p.
\end{align*}
We need an equation starting at $p$, not at $p g$, because membership in $\operatorname{Hol}_p(\omega)$ is defined using horizontal lifts that start at $p$. The principal equivariance of parallel transport supplies exactly this conversion:
\begin{align*}
\operatorname{PT}_{\gamma^{-1}}(p g)=\operatorname{PT}_{\gamma^{-1}}(p)g.
\end{align*}
Combining the two displayed equations gives
\begin{align*}
\operatorname{PT}_{\gamma^{-1}}(p)g=p.
\end{align*}
The right action of $G$ on the fiber $\pi^{-1}(x_0)$ is free and transitive. Therefore the unique element of $G$ carrying $\operatorname{PT}_{\gamma^{-1}}(p)$ to $p$ is $g$, equivalently
\begin{align*}
\operatorname{PT}_{\gamma^{-1}}(p)=p g^{-1}.
\end{align*}
Thus the reversed loop realizes $g^{-1}$ as a holonomy element based at $p$, so $g^{-1} \in \operatorname{Hol}_p(\omega)$.
[/guided]
[/step]
[step:Concatenate loops to obtain products]
Let $g_1,g_2 \in \operatorname{Hol}_p(\omega)$. Choose piecewise smooth loops
\begin{align*}
\gamma_1: [0,1] &\to M
\end{align*}
and
\begin{align*}
\gamma_2: [0,1] &\to M
\end{align*}
based at $x_0$ such that
\begin{align*}
\operatorname{PT}_{\gamma_1}(p)=p g_1
\end{align*}
and
\begin{align*}
\operatorname{PT}_{\gamma_2}(p)=p g_2.
\end{align*}
Let $\gamma_2 * \gamma_1: [0,1] \to M$ denote the piecewise smooth loop that first follows $\gamma_2$ and then follows $\gamma_1$. Parallel transport along a concatenated path is composition of the corresponding parallel transports in the same order as traversal, so
\begin{align*}
\operatorname{PT}_{\gamma_2 * \gamma_1}(p)=\operatorname{PT}_{\gamma_1}(\operatorname{PT}_{\gamma_2}(p)).
\end{align*}
Using the choice of $\gamma_2$ and then equivariance of $\operatorname{PT}_{\gamma_1}$,
\begin{align*}
\operatorname{PT}_{\gamma_2 * \gamma_1}(p)=\operatorname{PT}_{\gamma_1}(p g_2).
\end{align*}
\begin{align*}
\operatorname{PT}_{\gamma_1}(p g_2)=\operatorname{PT}_{\gamma_1}(p)g_2.
\end{align*}
Using the choice of $\gamma_1$ gives
\begin{align*}
\operatorname{PT}_{\gamma_2 * \gamma_1}(p)=(p g_1)g_2=p(g_1 g_2).
\end{align*}
Thus $g_1 g_2 \in \operatorname{Hol}_p(\omega)$.
[/step]
[step:Apply the subgroup criterion]
We have shown that $\operatorname{Hol}_p(\omega)$ contains the identity element $e$, is closed under inverses, and is closed under products. Therefore $\operatorname{Hol}_p(\omega)$ is a subgroup of $G$, that is,
\begin{align*}
\operatorname{Hol}_p(\omega) \le G.
\end{align*}
[/step]
