Advanced quantum technologies require scalable and controllable quantum resources1,2. Gaussian states of multimode light, such as squeezed states and cluster states, are scalable quantum systems3,4,5, which can be generated on demand. However, non-Gaussian features are indispensable in many quantum protocols, especially to reach a quantum computational advantage6. Embodying non-Gaussianity in a multimode quantum state remains a challenge as non-Gaussian operations generally cannot maintain coherence among multiple modes. Here, we generate non-Gaussian quantum states of a multimode light field by removing a single photon in a mode-selective manner from a Gaussian state7. To highlight the potential for continuous-variable quantum technologies, we first demonstrated the capability to generate negativity of the Wigner function in a controlled mode. Subsequently, we explored the interplay between non-Gaussianity and quantum entanglement and verify a theoretical prediction8 about the propagation of non-Gaussianity along the nodes of photon-subtracted cluster states. Our results demonstrate large-scale non-Gaussianity with great flexibility along with an ensured compatibility with quantum information protocols. This range of features makes our approach ideal to explore the physics of non-Gaussian entanglement9,10 and to develop quantum protocols, which range across quantum computing11,12, entanglement distillation13 and quantum simulations14.