Consider the finite dimensional unitary representations $\alpha,\beta,\gamma$ of the given compact group $G$ on corresponding vector spaces $V_1,V_2,V_3$. Let $|i\rangle_j,i=1,\dots,n_j$ be an orthnormal basis of $V_j$ where $dim V_j=n_j$. Then
$\{|i\rangle_1\otimes|j\rangle_2\otimes|k\rangle_3\}$ forms an orthonormal basis of $V=V_1\otimes V_2 \otimes V_3$. An element $g\in G$ acts on the tensor product space $V$ as
$$|i\rangle_1\otimes|j\rangle_2\otimes|k\rangle_3 \to \alpha(g) |i\rangle_1\otimes\beta(g)|j\rangle_2\otimes \gamma(g)|k\rangle_3 \tag 1$$
We can also find an orthogonal basis $e_{\mu},\mu=1,\dots,N$ (where $N=n_1n_2n_3$) of $V$ relative to which all elements $g \in G$ act as block diagonal matrices. More precisely, suppose with respect to basis $\{e_{\mu}\}$, the action of an element $g\in G$ on $V$ be denoted as
$$e_{\mu}\to \rho (g)e_\mu \tag 2$$
then $\rho(g)^{\nu}_{\mu}=\langle \nu|\rho(g)|\mu\rangle$ is a block diagonal matrix where the dimensions of different blocks$^1$ are independent of $g$. Let
$$|i\rangle_1\otimes|j\rangle_2\otimes|k\rangle_3=\sum_{\mu}\epsilon ^{\mu}_{ijk}e_{\mu}\tag 3$$
Acting with $g\in G$ on both sides of this equation gives
$$\alpha (g)|i\rangle_1\otimes \beta (g)|j\rangle_2\otimes \gamma(g)|k\rangle_3=\sum_{\mu}\epsilon ^{\mu}_{ijk}\rho (g)e_{\mu} \tag 4$$
Taking the scalar product with $|i'\rangle_1\otimes |j'\rangle_2 \otimes |k'\rangle_3$ and using (3) we get
$$\alpha(g)^{i'}_{i}\beta (g)^{j'}_{j}\gamma(g)^{k'}_{k}=\sum _{\mu,\nu} \rho^{\nu}_{\mu}(g)\epsilon^{*\nu}_{i'j'k'}{\epsilon}^{\mu}_{ijk}\tag 5$$
Now, according to the Peter-Weyl theorem (part 2), matrix elements of the irreducible representations of $G$ form an orthogonal basis of the space of square integrable functions on $G$ wrt the inner product
$$(A,B)=\int_G dg\; A(g)^{*}B(g) \tag 6$$
where $dg$ is the Haar measure. So, if we integrate both sides of (5), the nonzero contribution on RHS will only come from the part of $\rho$ which is direct sum of identity representations. In other words, let $W\subseteq V$ be the subspace of $V$ on which $G$ acts trivially, and let $\{e_1,\dots e_m\}\subseteq $ $\{e_1,\dots e_m,\dots,e_N\}$ be the basis of $W$, then the integration of (5) gives
$$\int_G dg\;\alpha(g)^{i'}_{i}\beta (g)^{j'}_{j}\gamma(g)^{k'}_{k}=\sum _{\mu,\nu =1}^{m} \delta^{\nu}_{\mu}\epsilon^{*\nu}_{i'j'k'}{\epsilon}^{\mu}_{ijk}=\sum _{\mu =1}^{m}\epsilon^{*\mu}_{i'j'k'}{\epsilon}^{\mu}_{ijk} \tag 7$$
where we have assumed that $Vol(G)=\displaystyle\int_G\; dg =1$
For the second part of your question, I would recommend these lecture notes. The basic idea for computing Wilson loop averages is following -
For a surface with boundary, the partition function of two dimensional Yang-Mills theory depends on the holonomy along the boundary. Let the partition function of a surface of genus $h$ and one boundary be denoted as $Z_h(U,ag^2)$ where $U$ is the fixed holonomy along the given boundary, $a$ is the area of the surface and $g$ is the Yang-Mills coupling constant. Now, consider the simplest situation in which a Wilson loop $W$ in representation $R_W$ is inserted along a contractible loop $C$ on a closed surface of genus $h$, and area $a$. To compute the Wilson loop average, we first cut the surface along $C$, which gives a disc $D$ of area (say) $b$ and another surface $S$ of area $c=a-b$, genus $h$ and one boundary. Now the Wilson loop average is given by integrating over $G$ the product of i) the partition functions of $D$ ii) partition function of $S$ and iii) the trace of the Wilson loop in representation $R_W$ -
$$\langle W\rangle=\frac {1}{Z_h(ag^2)}\int \: dU \:Z_h(U,(a-b)g^2)\chi_{R_W}(U)Z_0(U^{-1},bg^2) \tag 8$$
The case of a self-intersecting Wilson loop too is not very different.
$^1$ The smallest blocks form irreducible representations of $G$; Exactly which irreducible representations show up will depend on $\alpha,\beta,\gamma$; The same irreducible representations may also appear more than once
