Do you really want separate LSTM networks for each feature? Why not just have a single LSTM that takes a vector containing (X_1,t, X_2,t, …X_n,t) as input at time t and predicts a vector containing (X_1,t+1, X_2,t+1, …X_n,t+1, target_t+1)?
One advantage of this method is that the LSTM will be able to use the values all of the features when predicting the next value of each one, rather than only being able to use X_1,t<=t when predicting X_1,t+1.
If you want to do it your way, then you can make an LSTM model per feature that predicts the next feature value (and I would suggest that each model should also try to predict target_t+1), then the outputs of all the LSTM models could be concatenated and fed into a Linear layer or a simple feed-forward network in order to produce the combined prediction for target_t+1.
For the loss you just add them up as you suggest and it works as expected.
global_loss = loss1 + loss2 + ...
global_loss.backward()