We give an almost purely combinatorial expression for Wilson loop expectations of the Yang-Mills holonomy process with values in the unitary group on a compact oriented surface, possibly with boundary and arbitrary boundary conditions. Our main result computes the non-normalized expectation of products of traces of holonomies along an arbitrary family of immersed curves with transverse self-intersections and no triple points. It is expressed as a sum over assignments of highest weights of the unitary group to the connected components of the complement of the curves. Each term is a product of a Gaussian exponential factor, dimensions of unitary representations, and local contributions at the intersection points given by the sine or cosine of an angle determined by the surrounding highest weights. As an application, we obtain a new and short proof of the Makeenko-Migdal equations on arbitrary compact surfaces.