Pair microstructure of concentrated Brownian suspensions in simple-shear flow is studied by sampling of configurations from dynamic simulations by the Stokesian Dynamics technique. Simulated motions are three dimensional with periodic boundary conditions to mimic an infinitely extended suspension. Hydrodynamic interactions through Newtonian fluid and Brownian motion are the only physical influences upon the motion of the monodisperse hard-sphere particles. The dimensionless parameters characterizing the suspension are the particle volume fraction and Péclet number, defined, respectively, as φ=(4π/3)na3 with n the number density and a the sphere radius, and Pe=6πηγ̇a3/kT with η the fluid viscosity, γ̇ the shear rate, and kT the thermal energy. The majority of the results reported are from simulations at Pe=1000; results of simulations at Pe=1, 25, and 100 are also reported for φ=0.3 and φ=0.45. The pair structure is characterized by the pair distribution function, g(r)=P1|1(r)/n, where P1|1(r) is the conditional probability of finding a pair at a separation vector r. The structure under strong shearing exhibits an accumulation of pair probability at contact, and angular distortion (from spherical symmetry at Pe=0), with both effects increasing with Pe. Flow simulations were performed at Pe=1000 for eight volume fractions in the range 0.2⩽φ⩽0.585. For φ=0.2–0.3, the pair structure at contact, g(|r|=2)≡g(2), is found to exhibit a single region of strong correlation, g(2)≫1, at points around the axis of compression, with a particle-deficient wake in the extensional zones. A qualitative change in microstructure is observed between φ=0.3 and φ=0.37. For φ⩾0.37, the maximum g(2) lies at points in the shear plane nearly on the x axis of the bulk simple shear flow Ux=γ̇y, while at smaller φ, the maximum g(2) lies near the compressional axis; long-range string ordering is not observed. For φ=0.3 and φ=0.45, g(2)∼Pe0.7 for 1⩽Pe⩽1000, a slower increase than the g(2)∼Pe predicted theoretically for φ≪1 [Brady and Morris, J. Fluid Mech. 348, 143 (1997)]. Spherical harmonic decomposition of g(r) was performed, and for Pe=1000, expansion convergence is found to be nearly complete when harmonics Ylm to the level l=10 are included.