Abstract: In this paper, the immersed boundary method was used to simulate the vortex-induced vibration of a slender flexible cylinder exposed to linear shear flows. The vibration of the cylinder was simulated by using a three-dimensional cable model pinned at both ends. The cylinder has a mass ratio of 6 and an aspect ratio of 50. The normalized top tension is 496. The incoming flow has different linear velocity profiles with the shear rates ranging from 0 to 0.024. The maximum Reynold number is 250. It was found that the transverse vibration shows a standing wave pattern while the streamwise vibration shows a combined traveling-standing wave pattern. With the increase of the shear rate, the distribution of Power Spectrum Density (PSD) of vibration responses shows a multi-frequency mode, and the vibration energy shifts to low frequency bands. The mean drag coefficient varies in the spanwise direction while the root-mean-square (RMS) values of fluctuating drag and lift coefficients show a two-peak pattern. The distribution of the fluid-solid energy transferring coefficient indicates that the vibration-exciting region coincides with the high-velocity region while the vibration-damping region matches the low-velocity region. For the cases with low shear rates, the vortex-shedding behind the cylinder shows the interwoven pattern. However, for the cases with large shear rates, the vortex-shedding displays the oblique pattern. Due to the spanwise variation of the vortex-shedding frequency, the vortex-splitting occurs in the near-wake, leading to the vortex cells with different vortex-shedding frequencies.