  We consider a non-integrable model for interacting dimers on the
two-dimensional square lattice. Configurations are perfect matchings of
$\mathbb Z^2$, i.e. subsets of edges such that each vertex is covered exactly
once ("close-packing" condition). Dimer configurations are in bijection with
discrete height functions, defined on faces $\boldsymbol{\xi}$ of $\mathbb
Z^2$. The non-interacting model is "integrable" and solvable via Kasteleyn
theory; it is known that all the moments of the height difference
$h_{\boldsymbol{\xi}}-h_{\boldsymbol{\eta}}$ converge to those of the massless
Gaussian Free Field (GFF), asymptotically as
$|{\boldsymbol{\xi}}-{\boldsymbol{\eta}}|\to \infty$. We prove that the same
holds for small non-zero interactions, as was conjectured in the theoretical
physics literature. Remarkably, dimer-dimer correlation functions are instead
not universal and decay with a critical exponent that depends on the
interaction strength. Our proof is based on an exact representation of the
model in terms of lattice interacting fermions, which are studied by
constructive field theory methods. In the fermionic language, the height
difference $h_{\boldsymbol{\xi}}-h_{\boldsymbol{\eta}}$ takes the form of a
non-local operator, consisting of a sum of monomials along an {\it arbitrary}
path connecting $\boldsymbol{\xi}$ and $\boldsymbol{\eta}$. As in the
non-interacting case, this path-independence plays a crucial role in the proof.
