If very massive stars (M 100 M ⊙ ) can form and avoid too strong mass loss during their evolution, they are predicted to explode as pair-instability supernovae (PISNe). One critical test for candidate events is whether their nucleosynthesis yields and internal ejecta structure, being revealed through nebular-phase spectra at t 1 yr, match those of model predictions. Here we compute theoretical spectra based on model PISN ejecta at 1-3 years post-explosion to allow quantitative comparison with observations. The high column densities of PISNe lead to complete gamma-ray trapping for t 2 years which, combined with fulfilled conditions of steady state, leads to bolometric supernova luminosities matching the 56 Co decay. Most of the gamma-rays are absorbed by the deep-lying iron and silicon/sulphur layers. The ionization balance shows a predominantly neutral gas state, which leads to emission lines of Fe I, Si I, and S I. For low-mass PISNe the metal core expands slowly enough to produce a forest of distinct lines, whereas high-mass PISNe expand faster and produce more featureless spectra. Line blocking is complete below ∼5000Å for several years, and the model spectra are red. The strongest line is typically [Ca II] λλ7291, 7323, one of few lines from ionized species. We compare our models with proposed PISN candidates SN 2007bi and PTF12dam, finding discrepancies for several key observables and thus no support for a PISN interpretation. We discuss distinct spectral features predicted by the models, and the possibility of detecting pair-instability explosions among non-superluminous supernovae.