This work presents a novel nonlinear topology optimization formulation for designing soft pneumatic actuators. It combines density based topology optimization with porohyperelasticity in order to account, during the optimization, for the typically very large deformations occurring upon actuation. The pneumatic actuation is imposed as a predefined pressure source, applied to a fixed region of the design domain. In the entire design domain, the permeability field is defined as a function of the density field, approaching a large value in void regions, with density close to zero, and a small value in solid regions, with density close to one. The stiffness field is also defined as a function of the density field, using RAMP interpolation in order to penalize intermediate values of the density field. Portion of the boundary of the design domain is fixed, while an output region of the design domain is coupled to an external spring with a given stiffness, which represents the environment that the actuator should interact with. The displacement of the output region is maximized during the optimization. At the same time, a strain energy density constraint is imposed in order to avoid structures exhibiting extremely large strains under the applied pressure load. Unlike corresponding linear poroelasticity models, which are valid only for infinitesimal deformations, the present formulation produces different optimized actuator designs depending on the compliance level of the environment. Numerical examples demonstrate the capabilities of the proposed formulation in generating 3D designs which are manufacturable with additive manufacturing, and exhibit significantly superior performance compared to corresponding conventional designs. The efficient, accurate, and generic computational framework for soft pneumatic actuators, presented in this work, is envisioned to serve as a universal tool for designing devices such as gripper fingers and dexterous endoscopy tubes.