Another view which results in superior interpolation filter design is revealed by reviewing the spectrum of the signal with zeros inserted (resulting in an increase of the sampling rate by $I$ when you insert $I-1$ zeroes) and you will see that the original spectrum is replicated at integer multiples of the original sampling rate. The ideal low pass filter is the one that can pass the original spectrum with no distortion and reject all the new alias copies perfectly. (Ideal is not achievable but this defines the filter design targets to approach that). The ideal filter will interpolate the non zero values to grow the zeros to the value for the final interpolated waveform.
What is occurring, to explain the OP's question conceptually, is that the signal requires high frequency content when going from zeros to a large value in just one sample (and inserting zeros causes this to occur as evidenced by the additional high frequency content in the spectrum). This makes complete sense since the frequency content would be related to a change in magnitude versus a change in time (a big change in a short time requires high frequencies). A filter that removes the high frequencies results in the inability for the signal to change rapidly from one sample to the next (due to the memory of the previous samples which is what a filter does). A very easy way to see this occurring is with a simple moving average filter over $M$ samples when you interpolate by $M$ by inserting $M-1$ zeros; such a filter will hold the last sample as a zero-order hold, growing each zero to the last non-zero sample. This isn't a good way to do interpolation since only one non-zero sample is in the memory of the filter, but very easy to visualize. The "magic" happens with better designed filters as described in this answer which then consider many more samples to perform what is essentially interpolation with higher order polynomials.
With regards to Sinc reconstruction that is commonly used to explain "perfect" reconstruction; it is worth mentioning that Sinc reconstruction is only perfect when you can use a Sinc, but unfortunately the Sinc function extends to infinity which is impossible to do in practical realization (it is equivalent to filtering with a brickwall filter). Therefore the Sinc response is truncated in practice, which on its own results in a very poor filter, which is then typically windowed to improve performance. Ultimately, its use when truncated and windowed results in inferior reconstruction filters. Given the approach of paying attention to where the aliases actually reside in frequency, you can realize optimized filters (superior to simply truncating and windowing a Sinc function) using the multi-band filter algorithms such as least squares using
firls() function that is available in MATLAB, Python and Octave. I demonstrate this in the plots below from an X4 interpolation example with $f_s = 10$ KHz. The upper plot compares options for interpolation filters with the same number of taps, showing the poor performance of a truncated Sinc in green, and then the improved windowed Sinc filter in green along with the least squares multiband filter in red. All have the same number of taps, but the multiband approach results in 10 dB better rejection of the higher frequency images (and not visible from the scale of the plot, but lower passband ripple distortion as well).
The filter magnitude for the multiband filter was shifted down to overlay with the spectrum after zero-insert showing how the multiband filter maximizes rejection where needed thus resulting in the least amount of distortion in the resulting interpolated signal.
So the spectrum above with it's replication at the higher frequencies (images) has the zeros inserted in the time domain. When those higher frequencies are removed through filtering, the spectrum represents the original signal, at the higher sampling rate, as if it was sampled directly.
A truncated and properly windowed Sinc filter will not achieve the rejection or passband performance of this filter with the same number of taps since it would provide rejection at all locations in between the aliases, while the multi-band filter maximizes the rejection where it is needed. The resulting distortion is predictable and can be traded with the desired filter complexity. This is a standard interpolation design approach for high performance and efficient interpolation (especially when converted to polyphase structures!).