Is there any (open source) FEM toolbox that allows the direct discretization of higher order PDEs without the need to split them up into systems of second order?
FEniCS has a biharmonic example that does uses a mixed continuous-discontinuous Galerkin formulation. Any package with $C^0$ elements that can compute second derivatives and can integrate DG jump terms can also use this approach.
PetIGA Supports isogeometric elements with arbitrary continuity, but currently only on one "patch".
GeoPDEs appears to be capable, but I haven't used it.
OOFEM supports multi-patch isogeometric analysis with T-splines.
Note that it is not always better to use methods with higher continuity than to reduce to second order form. The bandwidth and vertex separators for these higher continuity methods are typically much larger, sometimes leading to more expensive solves than if you had reduced to second order form. You have to consider inherent heterogeneity/nonsmoothness in your problem, practical orders of accuracy, dimension, solution algorithm, and problem scale to determine which methods to prefer. Comparison studies using appropriate algorithms in each case would be a welcome contribution to the literature.
Have a look at the FENICS project. I have not used it yet, but it seems one can do a lot with it, with respect to finite element methods.