We investigate the effects of fiber aspect ratio, roughness, flexibility, and volume fraction on the rheology of concentrated suspensions in a steady shear flow using direct numerical simulations. We model the fibers as inextensible continuous flexible slender bodies with the Euler-Bernoulli beam equation governing their dynamics suspended in an incompressible Newtonian fluid. The fiber dynamics and fluid flow coupling is achieved using the immersed boundary method. In addition, the fiber surface roughness might lead to interfiber contacts, resulting in normal and tangential forces between the fibers, which follow Coulomb's law of friction. The surface roughness is modeled as hemispherical pro-trusions on the fiber surfaces. Their deformation results in a normal load-dependent friction coefficient. Our simulations accurately predict the experimentally observed shear thinning in fiber suspensions. Furthermore, we find that the suspension viscosity eta increases with increasing the volume fraction, roughness, fiber rigidity, and aspect ratio. The increase in eta is the macroscopic manifestation of a similar increase in the microscopic contact contribution to the total stress with these parameters. In addition, we observe positive and negative first N1 and second N2 normal stress differences, respectively, with |N2| < |N1|, in agreement with previous experiments. Last, we propose a modified Maron-Pierce law to quantify the reduction in the jamming volume fraction by increasing the fiber aspect ratio and roughness. Our results and analysis establish the use of fiber surface tribology to tune the suspension flow behavior.