We present a 3D numerical investigation of the viscometric functions for suspensions in viscoelastic fluids under steady shear including comparisons with theoretical predictions in the dilute regime and experimental results for non-colloidal suspensions in a Boger fluid. The comparison with dilute suspension theoretical predictions resolves a discrepancy in the literature regarding the second normal stress difference calculated using the Landau Lifshitz averaging procedure versus using the ensemble averaging procedure. We also determine that the particle contribution to the stress in a dilute suspension in viscoelastic fluids shear-thickens in all viscometric material properties and we present the scaling of the viscometric functions with Weissenberg number (Wi). This shear-thickening behavior is fundamentally different from that shown in 2D numerical simulations which are attributed to elongation flow effects between particles. Comparisons with experimental results at volume fraction phi = 0.05, which is the lowest volume fraction available in existing experimental measurements of bulk viscometric functions, highlight the important role of long range hydrodynamic interactions between particles to fully describe the suspension rheology even at low volume fractions. We investigate the effect of hydrodynamic interactions in the flow, vorticity, and gradient directions and show that these interactions can greatly enhance or decrease the viscometric functions as well as change the shear rate dependent behavior of the suspension. (C) 2016 Elsevier B.V. All rights reserved.