Seasonal differencing a when the sample size is small when estimating a VAR When estimating a VAR the series must be level and seasonally stationary. However I have only 48 data points. I first made the series level stationary based on the ADF test and performed HEGY test for seasonal unit roots. Few series which were level stationary exhibited seasonal unit roots. I am unable to take the seasonal difference due to the small sample size. How to deal with this issue?
I have included a plot of transformed series (by season) submitting to HEGY test below:

 A: Are you sure the series have seasonal unit roots? Given seasonal integration, neighboring observations (that obviously belong to different seasons) may diverge from each other as the sample size grows. Here is a simulated example of a seasonally-integrated quarterly series:

Is this what you are observing?
The null hypothesis of HEGY is presence of a seasonal unit root. Given a small sample, the power of the test is likely low, so it might just be struggling to reject a false $H_0$. Since your sample is only 4 observations per season, I would not trust the result of the test. While the test statistic is likely derived using an asymptotic approximation, 4 is quite far away from $\infty$... And this is not specific to HEGY alone; every test needs data to be powerful. (Well, there is an exception to every rule, but you get the general idea.)
I am not saying your series cannot have seasonal unit roots, but there is likely too little data to conclude anything much about it.

Regarding

When estimating a VAR the series must be level and seasonally stationary

see Ashley & Verbrugge "To difference or not to difference: a Monte Carlo investigation of inference in vector autoregression models" (2009). Here are some notes of the paper I took for myself:

*

*For model specification (short-term Granger causality) testing, use either a VAR in differences model or a lag-augmented VAR model (a levels VAR with lagged differences in addition to lagged levels); either has reasonably accurate size and comparable power in our simulations.

*Estimate the IRFs (and confidence intervals for them) using  a model in which the dependent variable is in levels – using either the levels model or the lag-augmented VAR model – with a trend term included.
The bias-corrected bootstrap confidence intervals appear to be somewhat preferable to using asymptotic standard errors. The actual coverage of nominally 95% confidence intervals may be substantially less than 95% except for very large samples, especially past lags one or two.

