In this paper, we establish a novel fractional model arising in the chemical reaction and develop an efficient spectral method for the three-dimensional Riesz-like space fractional nonlinear coupled reaction-diffusion equations. Based on the backward difference method for time stepping and the Legendre-Galerkin spectral method for space discretization, we construct a fully discrete numerical scheme which leads to a linear algebraic system. Then a direct method based on the matrix diagonalization approach is proposed to solve the linear algebraic system, where the cost of the algorithm is of a small multiple of $N^4$ ($N$ is the polynomial degree in each spatial coordinate) flops for each time level. In addition, the stability and convergence analysis are rigorously established. We obtain the optimal error estimate in space, and the results also show that the fully discrete scheme is unconditionally stable and convergent of order one in time. Furthermore, numerical experiments are presented to confirm the theoretical claims. As the applications of the proposed method, the fractional Gray-Scott model is solved to capture the pattern formation with an analysis of the properties of the fractional powers.