  The small-$x$ deep inelastic scattering in the saturation region is governed
by the non-linear evolution of Wilson-line operators. In the leading
logarithmic approximation it is given by the BK equation for the evolution of
color dipoles. In the next-to-leading order the BK equation gets contributions
from quark and gluon loops as well as from the tree gluon diagrams with
quadratic and cubic nonlinearities. We calculate the gluon contribution to
small-x evolution of Wilson lines (the quark part was obtained earlier).
