Biological membranes can be considered two-dimensional fluids with suspended integral membrane proteins (IMPs). We have calculated the effect of hydrodynamic interactions on the various diffusion coefficients of IMPs in lipid bilayers. The IMPs are modelled as hard cylinders of radius a immersed in a thin sheet of viscosity mu and thickness h bounded by a fluid of low viscosity mu'. We have ensemble averaged the N-body Stokes equations to the pair level and have renormalized them following the methods of Batchelor (1972) and Hinch (1977). The lengthscale for the hydrodynamic interactions is lambda a = mu h/mu', which is O(100a), and the slow decay of the interactions introduces new features in the renormalizations compared to the analogous analyses for three-dimensional suspensions of spheres. We have calculated the asymptotic limits for the short- and long-time tracer diffusivities, D-s and D-l, respectively, and for the gradient diffusivity, D-g, for phi much less than 1 and lambda much greater than 1, where phi is the IMP area fraction and lambda = mu h/(mu'a). The diffusivities are D-s/D-0 = 1 - 2 phi[1-(1+In(2)-9/32)/(In(lambda)-gamma)], D-l/D-0 = D-3/D-0-0.07/(In(lambda)-gamma), D-g/D-0 = 1+phi[-7+(6In(2)+7/16+0.37)/(In(lambda)-gamma], where D-0 is the diffusivity in the limit of zero area fraction, and gamma = 0.577216 is Euler's constant. The results for D-l and D-s differ only slightly. The decrease in D-g/D-0 as phi increases contrasts with the result for spheres for which D-g/D-0 > 1.