We present a fluctuating boundary integral method (FBIM) for Brownian Dynamics (BD) of suspensions of rigid particles of complex shape immersed in a Stokes fluid. Our approach relies on a first-kind boundary integral formulation of a Stochastic Stokes Boundary Value Problem (SSBVP) in which a random surface velocity is prescribed on the particle surface. This random surface velocity has zero mean and covariance proportional to the Green’s function for the Stokes flow (Stokeslet). Furthermore, we demonstrate that discretizing the first-kind formulation using standard boundary integral techniques leads to an efficient numerical method that strictly preserves discrete fluctuation-dissipation balance (DFDB). We develop fast linear-scaling algorithms for performing matrix-vector products and for rapidly generating the random surface velocity by employing the Hasimoto splitting of the Stokeslet, which guarantees that the near-field (short-ranged) and far-field (long-ranged) contributions of the Stokeslet are independently symmetric and positive-definite. This ensures that the computaitonal cost of Brownian simulation is only marginally larger than the cost of deterministic simulations, in stark contrast to traditional BD approaches. To handle the inherent ill-conditioning of the linear system due to the first-kind formulation, we employ an effective block-diagonal preconditioner that ignores all the hydrodynamic interactions between distinct particles. FBIM provides the key ingredient for time integration of the overdamped Langevin equations for Brownian suspensions of rigid particles. We demonstrate that FBIM obeys DFDB by performing equilibrium BD simulations of suspensions of starfish-shaped bodies using a random finite difference temporal integrator. We numerically demonstrate the linear scaling of FBIM and discuss the computational cost of the various components of the algorithm.