Large intelligent surface (LIS) has recently emerged as a potential enabling technology for 6G networks, offering extended coverage and enhanced energy and spectral efficiency. In this work, motivated by its promising potentials, we investigate the error rate performance of LIS-assisted non-orthogonal multiple access (NOMA) networks. Specifically, we consider a downlink NOMA system, in which data transmission between a base station (BS) and L NOMA users is assisted by an LIS comprising M reflective elements. First, we derive the probability density function of the end-to-end wireless fading channels between the BS and NOMA users. Then, by leveraging the obtained results, we derive an approximate expression for the pairwise error probability (PEP) of NOMA users under the assumption of imperfect successive interference cancellation. Furthermore, accurate expressions for the PEP for M = 1 and large M values (M > 10) are presented in closed-form. To gain further insights into the system performance, an asymptotic expression for the PEP in high signal-to-noise ratio regime, the asymptotic diversity order, and a tight union bound on the bit error rate are provided. Finally, numerical and simulation results are presented to validate the derived mathematical results. Index terms-Asymptotic diversity order, bit error rate, large intelligent surfaces, NOMA, pairwise error probability, union bound, 6G.