If each individual subsystem can have a range of energy level
I'm not sure what you mean here. The summation in the partition function applies when the energy states are discrete. When they are continuous, as is the case in most of classical stat mech$^*$, the sum is replaced by an integral. This has nothing to do with subsystems; it is merely a consequence of switching to a continuous energy spectrum.
As @alarge said, the Boltzmann factor should be weighted by the degeneracy in the discrete case or the density of states in the continuous case. Then your partition function, to within a constant factor, is
$$Z = \int \Omega\left( E\right)e^{-E/kT}\mathrm{d}E$$
However, in most cases the calculation of $\Omega\left( E\right)$ is non-trivial at best, and outright impossible at worst.
A major benefit of the canonical ensemble is that we don't actually need to calculate $\Omega\left( E\right)$ if we know how the energy depends on the degrees of freedom within the system. For this we use the Hamiltonian function, which for, e.g. a gas particle, describes the energy in terms of the position and momentum of the particle. It is then sufficient to integrate the Boltzmann factor - with the Hamiltonian replacing $E$ - over all of phase space (all possible positions and momenta), like so:
$$Z = \int e^{-H\left(q, p\right)/kT}dqdp$$
where $q$ and $p$ are the position and momentum, respectively, and constants have again been omitted for clarity.
So to answer your question, it is not sufficient to integrate the energy over all states. If you want to integrate the energy, you need to weight it with the DOS. However, I highly suggest you take the approach of summing or integrating the Hamiltonian over all degrees of freedom.
$^*$ Discrete systems are generally easier to deal with, so they are often introduced first.
